1 Setup

1.1 Libs

#install.packages("remotes")
#install.packages('TMB',source=TRUE) #didn't work on home laptop
#remotes::install_github("glmmTMB/glmmTMB/glmmTMB")
library("glmmTMB")
## Warning in check_dep_version(dep_pkg = "TMB"): package version mismatch: 
## glmmTMB was built with TMB package version 1.9.15
## Current TMB package version is 1.9.16
## Please re-install glmmTMB from source or restore original 'TMB' package (see '?reinstalling' for more information)
#install.packages("BiocManager")
#BiocManager::install("phyloseq")
library("phyloseq")
library("dplyr")
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
# #install.packages("reshape")
library("reshape")
## 
## Attaching package: 'reshape'
## The following object is masked from 'package:dplyr':
## 
##     rename
# library("ggplot2")
# #install.packages("bbmle")
library("bbmle")
## Loading required package: stats4
## 
## Attaching package: 'bbmle'
## The following object is masked from 'package:dplyr':
## 
##     slice
# #install.packages("car")
# library("car")
library("DHARMa")
## This is DHARMa 0.4.7. For overview type '?DHARMa'. For recent changes, type news(package = 'DHARMa')
library("plyr")
## ------------------------------------------------------------------------------
## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)
## ------------------------------------------------------------------------------
## 
## Attaching package: 'plyr'
## The following objects are masked from 'package:reshape':
## 
##     rename, round_any
## The following objects are masked from 'package:dplyr':
## 
##     arrange, count, desc, failwith, id, mutate, rename, summarise,
##     summarize
library("ggplot2")

1.2 Re-read data

setwd("~/nicolagk@hawaii.edu - Google Drive/My Drive/Bamboo_mesos/Bamboo_mesos/04.mixed_models/")

#ps.clean <- readRDS("../01.process_asvs/reusum23.ps.lulu.clean.rds")
ps.clean <- readRDS("../01.process_asvs/bamboomesos_ps.lulu.clean.decontam.rds")
#ps.clean <- readRDS("../01.process_asvs/ps.all.clean.100.rds")
ps.clean #514 taxa, 169 samples
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 514 taxa and 169 samples ]
## sample_data() Sample Data:       [ 169 samples by 24 sample variables ]
## tax_table()   Taxonomy Table:    [ 514 taxa by 10 taxonomic ranks ]
ps.exp <- subset_samples(ps.clean,Mesocosm_type=="Experiment"&Microbe_treatment!="AHK"&Exp_day!=20&Microbe_treatment!="ECO"&Microbe_treatment!="ALL")
ps.exp #60 samples
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 514 taxa and 60 samples ]
## sample_data() Sample Data:       [ 60 samples by 24 sample variables ]
## tax_table()   Taxonomy Table:    [ 514 taxa by 10 taxonomic ranks ]
##including ECO & ALL to compare
ps.exp.all <- subset_samples(ps.clean,Mesocosm_type=="Experiment"&Microbe_treatment!="AHK"&Exp_day!=20)
ps.exp.all #100 samples
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 514 taxa and 100 samples ]
## sample_data() Sample Data:       [ 100 samples by 24 sample variables ]
## tax_table()   Taxonomy Table:    [ 514 taxa by 10 taxonomic ranks ]
##water samples
ps.exp.wtr <- subset_samples(ps.exp,Sample_type=="Water")
ps.exp.wtr #30 samples
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 514 taxa and 30 samples ]
## sample_data() Sample Data:       [ 30 samples by 24 sample variables ]
## tax_table()   Taxonomy Table:    [ 514 taxa by 10 taxonomic ranks ]
##plus eco & all
ps.exp.all.wtr <- subset_samples(ps.exp.all,Sample_type=="Water")
ps.exp.all.wtr #50 samples
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 514 taxa and 50 samples ]
## sample_data() Sample Data:       [ 50 samples by 24 sample variables ]
## tax_table()   Taxonomy Table:    [ 514 taxa by 10 taxonomic ranks ]
##larvae samples
ps.exp.lar <- subset_samples(ps.exp,Sample_type=="Larvae")
ps.exp.lar #30 samples
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 514 taxa and 30 samples ]
## sample_data() Sample Data:       [ 30 samples by 24 sample variables ]
## tax_table()   Taxonomy Table:    [ 514 taxa by 10 taxonomic ranks ]
##plus eco & all
ps.exp.all.lar <- subset_samples(ps.exp.all,Sample_type=="Larvae")
ps.exp.all.lar #50 samples
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 514 taxa and 50 samples ]
## sample_data() Sample Data:       [ 50 samples by 24 sample variables ]
## tax_table()   Taxonomy Table:    [ 514 taxa by 10 taxonomic ranks ]

2 Data formatting

2.1 Trim, thoroughly

Ran once then saved

##6 samples out of 50 = 12%
##6 samples out of 30 = 20%

##water
ps.exp.wtr.trim <- filter_taxa(ps.exp.wtr, function(x) sum(x > 1) > (0.20*length(x)), TRUE)
ps.exp.wtr.trim ##62 taxa
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 62 taxa and 30 samples ]
## sample_data() Sample Data:       [ 30 samples by 24 sample variables ]
## tax_table()   Taxonomy Table:    [ 62 taxa by 10 taxonomic ranks ]
tail(sort(sample_sums(ps.exp.wtr.trim),decreasing=T))
##   WB9  WL13   WB7   WL6  WL23   WL9 
## 77861 72372 68692 63574 53741 51123
##with all and eco
ps.exp.all.wtr.trim <- filter_taxa(ps.exp.all.wtr, function(x) sum(x > 1) > (0.12*length(x)), TRUE)
ps.exp.all.wtr.trim #69 taxa
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 69 taxa and 50 samples ]
## sample_data() Sample Data:       [ 50 samples by 24 sample variables ]
## tax_table()   Taxonomy Table:    [ 69 taxa by 10 taxonomic ranks ]
##larvae
ps.exp.lar.trim <- filter_taxa(ps.exp.lar, function(x) sum(x > 1) > (0.20*length(x)), TRUE)
ps.exp.lar.trim #44 taxa 
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 44 taxa and 30 samples ]
## sample_data() Sample Data:       [ 30 samples by 24 sample variables ]
## tax_table()   Taxonomy Table:    [ 44 taxa by 10 taxonomic ranks ]
##with all and eco
ps.exp.all.lar.trim <- filter_taxa(ps.exp.all.lar, function(x) sum(x > 1) > (0.12*length(x)), TRUE)
ps.exp.all.lar.trim #54 taxa
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 54 taxa and 50 samples ]
## sample_data() Sample Data:       [ 50 samples by 24 sample variables ]
## tax_table()   Taxonomy Table:    [ 54 taxa by 10 taxonomic ranks ]
tail(sort(sample_sums(ps.exp.lar.trim),decreasing=T))
##   LB8   LB9  LL12   LL7  LB11  LB12 
## 36318 34226 33371 31971 18211 13437
##proportion of zeroes
##water
otu.exp.wtr <- data.frame(ps.exp.wtr@otu_table)
sum(otu.exp.wtr == 0) / (ncol(otu.exp.wtr) * nrow(otu.exp.wtr)) * 100
## [1] 93.19066
#93.19066%
 
otu.exp.wtr.trim <- data.frame(ps.exp.wtr.trim@otu_table)
sum(otu.exp.wtr.trim == 0) / (ncol(otu.exp.wtr.trim) * nrow(otu.exp.wtr.trim)) * 100
## [1] 54.94624
#54.94624%

##water - with all
otu.exp.all.wtr <- data.frame(ps.exp.all.wtr@otu_table)
sum(otu.exp.all.wtr == 0) / (ncol(otu.exp.all.wtr) * nrow(otu.exp.all.wtr)) * 100
## [1] 94.27237
#94.27237%
 
otu.exp.all.wtr.trim <- data.frame(ps.exp.all.wtr.trim@otu_table)
sum(otu.exp.all.wtr.trim == 0) / (ncol(otu.exp.all.wtr.trim) * nrow(otu.exp.all.wtr.trim)) * 100
## [1] 64.11594
#64.11594%

##save if ready
#saveRDS(ps.exp.wtr.trim,file="ps.exp.wtr.trim.rds")
#saveRDS(ps.exp.all.wtr.trim,file="ps.exp.all.wtr.trim.rds")

##larvae
otu.exp.lar <- data.frame(ps.exp.lar@otu_table)
sum(otu.exp.lar == 0) / (ncol(otu.exp.lar) * nrow(otu.exp.lar)) * 100
## [1] 95.32425
#95.32425%
 
otu.exp.lar.trim <- data.frame(ps.exp.lar.trim@otu_table)
sum(otu.exp.lar.trim == 0) / (ncol(otu.exp.lar.trim) * nrow(otu.exp.lar.trim)) * 100
## [1] 61.81818
#61.81818%

##larvae - with all
otu.exp.all.lar <- data.frame(ps.exp.all.lar@otu_table)
sum(otu.exp.all.lar == 0) / (ncol(otu.exp.all.lar) * nrow(otu.exp.all.lar)) * 100
## [1] 95.96498
#95.96498%
 
otu.exp.all.lar.trim <- data.frame(ps.exp.all.lar.trim@otu_table)
sum(otu.exp.all.lar.trim == 0) / (ncol(otu.exp.all.lar.trim) * nrow(otu.exp.all.lar.trim)) * 100
## [1] 70.96296
#70.96296%

##save if ready
#saveRDS(ps.exp.lar.trim,file="ps.exp.lar.trim.rds")
#saveRDS(ps.exp.all.lar.trim,file="ps.exp.all.lar.trim.rds")

2.1.1 Which ASVs to subset from fasta

Didn’t do yet

otus.lar.trim <- colnames(otu.exp.lar.trim)
otus.wtr.trim <- colnames(otu.exp.wtr.trim)

#write.table(otus.lar.trim, file="otus.larvae.trim.txt", append = FALSE, sep = "/n", row.names = FALSE, col.names = FALSE,quote=FALSE)
#write.table(otus.wtr.trim, file="otus.water.trim.txt", append = FALSE, sep = "/n", row.names = FALSE, col.names = FALSE,quote=FALSE)

2.2 More formatting - long form

2.2.1 Just EMO, MMO, and MEM

otu.wtr <- data.frame(ps.exp.wtr.trim@otu_table)
otu.lar <- data.frame(ps.exp.lar.trim@otu_table)

#total read count colums
otu.wtr$sum <- rowSums(otu.wtr)
otu.lar$sum <- rowSums(otu.lar)

otu.wtr$Short_label <- row.names(otu.wtr)
otu.lar$Short_label <- row.names(otu.lar)

#bring in the metadata again
samdf.wtr <- data.frame(ps.exp.wtr.trim@sam_data)
samdf.lar <- data.frame(ps.exp.lar.trim@sam_data)

##add extra metadata
#meso.data <- read.csv("../bamboomesos_metadata.csv")

#samdf.wtr.more <- merge(samdf.wtr,meso.data)
#samdf.lar.more <- merge(samdf.lar,meso.data)

#samdf.wtr$food.microbes <- paste0(samdf.wtr$Food_type,"_",samdf.wtr$Microbe_treatment)
#samdf.lar$food.microbes <- paste0(samdf.lar$Food_type,"_",samdf.lar$Microbe_treatment)
# 
# #join the data sets by sample name
df.water <- plyr::join(otu.wtr, samdf.wtr, by="Short_label", type="left")
df.larv <- plyr::join(otu.lar, samdf.lar, by="Short_label", type="left")
colnames(df.water)
##  [1] "ASV0002"            "ASV0003"            "ASV0004"           
##  [4] "ASV0005"            "ASV0007"            "ASV0008"           
##  [7] "ASV0009"            "ASV0010"            "ASV0011"           
## [10] "ASV0012"            "ASV0015"            "ASV0016"           
## [13] "ASV0017"            "ASV0020"            "ASV0021"           
## [16] "ASV0022"            "ASV0023"            "ASV0024"           
## [19] "ASV0025"            "ASV0027"            "ASV0036"           
## [22] "ASV0037"            "ASV0039"            "ASV0040"           
## [25] "ASV0042"            "ASV0044"            "ASV0051"           
## [28] "ASV0054"            "ASV0057"            "ASV0059"           
## [31] "ASV0063"            "ASV0065"            "ASV0066"           
## [34] "ASV0069"            "ASV0072"            "ASV0073"           
## [37] "ASV0079"            "ASV0080"            "ASV0082"           
## [40] "ASV0084"            "ASV0085"            "ASV0089"           
## [43] "ASV0091"            "ASV0093"            "ASV0094"           
## [46] "ASV0095"            "ASV0101"            "ASV0102"           
## [49] "ASV0105"            "ASV0111"            "ASV0112"           
## [52] "ASV0115"            "ASV0116"            "ASV0119"           
## [55] "ASV0123"            "ASV0128"            "ASV0139"           
## [58] "ASV0145"            "ASV0146"            "ASV0155"           
## [61] "ASV0167"            "ASV0180"            "sum"               
## [64] "Short_label"        "Longer_name"        "Sample_type"       
## [67] "Num_larvae"         "Date_collected"     "Exp_day"           
## [70] "Mesocosm_id"        "Mesocosm_treatment" "Microbe_treatment" 
## [73] "Food_type"          "Mesocosm_type"      "Treatment_notes"   
## [76] "Notes"              "Raw_reads"          "Num_larvae_d00"    
## [79] "Num_larvae_d08"     "Num_larvae_d34"     "Num_larvae_d40"    
## [82] "Total_pupae"        "Proportion_pupated" "Biomass_day30"     
## [85] "Biomass_day40"      "lib_size_clean"     "is.neg"
##getting rid of some irrelevant metadat
colnames(df.water)
##  [1] "ASV0002"            "ASV0003"            "ASV0004"           
##  [4] "ASV0005"            "ASV0007"            "ASV0008"           
##  [7] "ASV0009"            "ASV0010"            "ASV0011"           
## [10] "ASV0012"            "ASV0015"            "ASV0016"           
## [13] "ASV0017"            "ASV0020"            "ASV0021"           
## [16] "ASV0022"            "ASV0023"            "ASV0024"           
## [19] "ASV0025"            "ASV0027"            "ASV0036"           
## [22] "ASV0037"            "ASV0039"            "ASV0040"           
## [25] "ASV0042"            "ASV0044"            "ASV0051"           
## [28] "ASV0054"            "ASV0057"            "ASV0059"           
## [31] "ASV0063"            "ASV0065"            "ASV0066"           
## [34] "ASV0069"            "ASV0072"            "ASV0073"           
## [37] "ASV0079"            "ASV0080"            "ASV0082"           
## [40] "ASV0084"            "ASV0085"            "ASV0089"           
## [43] "ASV0091"            "ASV0093"            "ASV0094"           
## [46] "ASV0095"            "ASV0101"            "ASV0102"           
## [49] "ASV0105"            "ASV0111"            "ASV0112"           
## [52] "ASV0115"            "ASV0116"            "ASV0119"           
## [55] "ASV0123"            "ASV0128"            "ASV0139"           
## [58] "ASV0145"            "ASV0146"            "ASV0155"           
## [61] "ASV0167"            "ASV0180"            "sum"               
## [64] "Short_label"        "Longer_name"        "Sample_type"       
## [67] "Num_larvae"         "Date_collected"     "Exp_day"           
## [70] "Mesocosm_id"        "Mesocosm_treatment" "Microbe_treatment" 
## [73] "Food_type"          "Mesocosm_type"      "Treatment_notes"   
## [76] "Notes"              "Raw_reads"          "Num_larvae_d00"    
## [79] "Num_larvae_d08"     "Num_larvae_d34"     "Num_larvae_d40"    
## [82] "Total_pupae"        "Proportion_pupated" "Biomass_day30"     
## [85] "Biomass_day40"      "lib_size_clean"     "is.neg"
df.water1 <- df.water %>%
  dplyr::select(-Short_label,-Longer_name,-Sample_type,-Num_larvae,-Date_collected,-Exp_day,-Mesocosm_treatment,-Mesocosm_type,-Treatment_notes,-Notes,-Raw_reads,-Num_larvae_d00,-Num_larvae_d08,-Num_larvae_d34,-Num_larvae_d40,-Total_pupae,-Proportion_pupated,-Biomass_day30,-lib_size_clean,-is.neg)
colnames(df.water1)
##  [1] "ASV0002"           "ASV0003"           "ASV0004"          
##  [4] "ASV0005"           "ASV0007"           "ASV0008"          
##  [7] "ASV0009"           "ASV0010"           "ASV0011"          
## [10] "ASV0012"           "ASV0015"           "ASV0016"          
## [13] "ASV0017"           "ASV0020"           "ASV0021"          
## [16] "ASV0022"           "ASV0023"           "ASV0024"          
## [19] "ASV0025"           "ASV0027"           "ASV0036"          
## [22] "ASV0037"           "ASV0039"           "ASV0040"          
## [25] "ASV0042"           "ASV0044"           "ASV0051"          
## [28] "ASV0054"           "ASV0057"           "ASV0059"          
## [31] "ASV0063"           "ASV0065"           "ASV0066"          
## [34] "ASV0069"           "ASV0072"           "ASV0073"          
## [37] "ASV0079"           "ASV0080"           "ASV0082"          
## [40] "ASV0084"           "ASV0085"           "ASV0089"          
## [43] "ASV0091"           "ASV0093"           "ASV0094"          
## [46] "ASV0095"           "ASV0101"           "ASV0102"          
## [49] "ASV0105"           "ASV0111"           "ASV0112"          
## [52] "ASV0115"           "ASV0116"           "ASV0119"          
## [55] "ASV0123"           "ASV0128"           "ASV0139"          
## [58] "ASV0145"           "ASV0146"           "ASV0155"          
## [61] "ASV0167"           "ASV0180"           "sum"              
## [64] "Mesocosm_id"       "Microbe_treatment" "Food_type"        
## [67] "Biomass_day40"
df.water.long <- melt(df.water1, id.vars=c("sum", "Mesocosm_id", "Microbe_treatment","Food_type", "Biomass_day40"))

df.water.long$rowid <- 1:nrow(df.water.long)

df.larv1 <- df.larv %>%
  dplyr::select(-Short_label,-Longer_name,-Sample_type,-Num_larvae,-Date_collected,-Exp_day,-Mesocosm_treatment,-Mesocosm_type,-Treatment_notes,-Notes,-Raw_reads,-Num_larvae_d00,-Num_larvae_d08,-Num_larvae_d34,-Num_larvae_d40,-Total_pupae,-Proportion_pupated,-Biomass_day30,-lib_size_clean,-is.neg)
colnames(df.larv1)
##  [1] "ASV0002"           "ASV0003"           "ASV0005"          
##  [4] "ASV0006"           "ASV0007"           "ASV0009"          
##  [7] "ASV0010"           "ASV0011"           "ASV0012"          
## [10] "ASV0013"           "ASV0016"           "ASV0019"          
## [13] "ASV0024"           "ASV0027"           "ASV0032"          
## [16] "ASV0033"           "ASV0043"           "ASV0053"          
## [19] "ASV0054"           "ASV0055"           "ASV0058"          
## [22] "ASV0059"           "ASV0060"           "ASV0061"          
## [25] "ASV0066"           "ASV0070"           "ASV0079"          
## [28] "ASV0083"           "ASV0085"           "ASV0087"          
## [31] "ASV0099"           "ASV0100"           "ASV0102"          
## [34] "ASV0125"           "ASV0130"           "ASV0134"          
## [37] "ASV0138"           "ASV0149"           "ASV0150"          
## [40] "ASV0164"           "ASV0218"           "ASV0262"          
## [43] "ASV0289"           "ASV0311"           "sum"              
## [46] "Mesocosm_id"       "Microbe_treatment" "Food_type"        
## [49] "Biomass_day40"
#df.larv2 <- df.larv1[df.larv1$sum!=0,]
#df.larv2$sum==0

df.larv.long <- melt(df.larv1, id.vars=c("sum", "Mesocosm_id", "Microbe_treatment","Food_type", "Biomass_day40"))

df.larv.long$rowid <- 1:nrow(df.larv.long)

##save
#saveRDS(df.water.long,file="df.water.long.rds")
#saveRDS(df.larv.long,file="df.larv.long.rds")

2.2.2 All of ’em

otu.all.wtr <- data.frame(ps.exp.all.wtr.trim@otu_table)
otu.all.lar <- data.frame(ps.exp.all.lar.trim@otu_table)

#total read count colums
otu.all.wtr$sum <- rowSums(otu.all.wtr)
otu.all.lar$sum <- rowSums(otu.all.lar)

otu.all.wtr$Short_label <- row.names(otu.all.wtr)
otu.all.lar$Short_label <- row.names(otu.all.lar)

#bring in the metadata again
samdf.all.wtr <- data.frame(ps.exp.all.wtr.trim@sam_data)
samdf.all.lar <- data.frame(ps.exp.all.lar.trim@sam_data)

##add extra metadata
#meso.data <- read.csv("../bamboomesos_metadata.csv")

#samdf.wtr.more <- merge(samdf.wtr,meso.data)
#samdf.lar.more <- merge(samdf.lar,meso.data)

#samdf.wtr$food.microbes <- paste0(samdf.wtr$Food_type,"_",samdf.wtr$Microbe_treatment)
#samdf.lar$food.microbes <- paste0(samdf.lar$Food_type,"_",samdf.lar$Microbe_treatment)
# 
# #join the data sets by sample name
df.all.water <- plyr::join(otu.all.wtr, samdf.all.wtr, by="Short_label", type="left")
df.all.larv <- plyr::join(otu.all.lar, samdf.all.lar, by="Short_label", type="left")
colnames(df.all.water)
##  [1] "ASV0001"            "ASV0002"            "ASV0003"           
##  [4] "ASV0004"            "ASV0005"            "ASV0007"           
##  [7] "ASV0008"            "ASV0009"            "ASV0010"           
## [10] "ASV0011"            "ASV0012"            "ASV0015"           
## [13] "ASV0016"            "ASV0017"            "ASV0018"           
## [16] "ASV0020"            "ASV0021"            "ASV0022"           
## [19] "ASV0023"            "ASV0024"            "ASV0025"           
## [22] "ASV0027"            "ASV0033"            "ASV0036"           
## [25] "ASV0037"            "ASV0039"            "ASV0040"           
## [28] "ASV0042"            "ASV0044"            "ASV0051"           
## [31] "ASV0054"            "ASV0057"            "ASV0059"           
## [34] "ASV0063"            "ASV0065"            "ASV0066"           
## [37] "ASV0069"            "ASV0072"            "ASV0073"           
## [40] "ASV0079"            "ASV0080"            "ASV0082"           
## [43] "ASV0083"            "ASV0084"            "ASV0085"           
## [46] "ASV0089"            "ASV0091"            "ASV0093"           
## [49] "ASV0094"            "ASV0095"            "ASV0101"           
## [52] "ASV0102"            "ASV0105"            "ASV0111"           
## [55] "ASV0112"            "ASV0115"            "ASV0116"           
## [58] "ASV0119"            "ASV0123"            "ASV0128"           
## [61] "ASV0139"            "ASV0145"            "ASV0146"           
## [64] "ASV0152"            "ASV0155"            "ASV0167"           
## [67] "ASV0175"            "ASV0180"            "ASV0184"           
## [70] "sum"                "Short_label"        "Longer_name"       
## [73] "Sample_type"        "Num_larvae"         "Date_collected"    
## [76] "Exp_day"            "Mesocosm_id"        "Mesocosm_treatment"
## [79] "Microbe_treatment"  "Food_type"          "Mesocosm_type"     
## [82] "Treatment_notes"    "Notes"              "Raw_reads"         
## [85] "Num_larvae_d00"     "Num_larvae_d08"     "Num_larvae_d34"    
## [88] "Num_larvae_d40"     "Total_pupae"        "Proportion_pupated"
## [91] "Biomass_day30"      "Biomass_day40"      "lib_size_clean"    
## [94] "is.neg"
##getting rid of some irrelevant metadat
df.all.water1 <- df.all.water %>%
  dplyr::select(-Short_label,-Longer_name,-Sample_type,-Num_larvae,-Date_collected,-Exp_day,-Mesocosm_treatment,-Mesocosm_type,-Treatment_notes,-Notes,-Raw_reads,-Num_larvae_d00,-Num_larvae_d08,-Num_larvae_d34,-Num_larvae_d40,-Total_pupae,-Proportion_pupated,-Biomass_day30,-lib_size_clean,-is.neg)
colnames(df.all.water1)
##  [1] "ASV0001"           "ASV0002"           "ASV0003"          
##  [4] "ASV0004"           "ASV0005"           "ASV0007"          
##  [7] "ASV0008"           "ASV0009"           "ASV0010"          
## [10] "ASV0011"           "ASV0012"           "ASV0015"          
## [13] "ASV0016"           "ASV0017"           "ASV0018"          
## [16] "ASV0020"           "ASV0021"           "ASV0022"          
## [19] "ASV0023"           "ASV0024"           "ASV0025"          
## [22] "ASV0027"           "ASV0033"           "ASV0036"          
## [25] "ASV0037"           "ASV0039"           "ASV0040"          
## [28] "ASV0042"           "ASV0044"           "ASV0051"          
## [31] "ASV0054"           "ASV0057"           "ASV0059"          
## [34] "ASV0063"           "ASV0065"           "ASV0066"          
## [37] "ASV0069"           "ASV0072"           "ASV0073"          
## [40] "ASV0079"           "ASV0080"           "ASV0082"          
## [43] "ASV0083"           "ASV0084"           "ASV0085"          
## [46] "ASV0089"           "ASV0091"           "ASV0093"          
## [49] "ASV0094"           "ASV0095"           "ASV0101"          
## [52] "ASV0102"           "ASV0105"           "ASV0111"          
## [55] "ASV0112"           "ASV0115"           "ASV0116"          
## [58] "ASV0119"           "ASV0123"           "ASV0128"          
## [61] "ASV0139"           "ASV0145"           "ASV0146"          
## [64] "ASV0152"           "ASV0155"           "ASV0167"          
## [67] "ASV0175"           "ASV0180"           "ASV0184"          
## [70] "sum"               "Mesocosm_id"       "Microbe_treatment"
## [73] "Food_type"         "Biomass_day40"
df.all.water.long <- melt(df.all.water1, id.vars=c("sum", "Mesocosm_id", "Microbe_treatment","Food_type", "Biomass_day40"))

df.all.water.long$rowid <- 1:nrow(df.all.water.long)

df.all.larv1 <- df.all.larv %>%
  dplyr::select(-Short_label,-Longer_name,-Sample_type,-Num_larvae,-Date_collected,-Exp_day,-Mesocosm_treatment,-Mesocosm_type,-Treatment_notes,-Notes,-Raw_reads,-Num_larvae_d00,-Num_larvae_d08,-Num_larvae_d34,-Num_larvae_d40,-Total_pupae,-Proportion_pupated,-Biomass_day30,-lib_size_clean,-is.neg)
colnames(df.all.larv1)
##  [1] "ASV0001"           "ASV0002"           "ASV0003"          
##  [4] "ASV0004"           "ASV0005"           "ASV0006"          
##  [7] "ASV0007"           "ASV0009"           "ASV0010"          
## [10] "ASV0011"           "ASV0012"           "ASV0013"          
## [13] "ASV0016"           "ASV0019"           "ASV0024"          
## [16] "ASV0027"           "ASV0030"           "ASV0032"          
## [19] "ASV0033"           "ASV0043"           "ASV0053"          
## [22] "ASV0054"           "ASV0055"           "ASV0058"          
## [25] "ASV0059"           "ASV0060"           "ASV0061"          
## [28] "ASV0066"           "ASV0070"           "ASV0079"          
## [31] "ASV0083"           "ASV0085"           "ASV0087"          
## [34] "ASV0099"           "ASV0100"           "ASV0102"          
## [37] "ASV0109"           "ASV0125"           "ASV0130"          
## [40] "ASV0134"           "ASV0138"           "ASV0149"          
## [43] "ASV0150"           "ASV0152"           "ASV0164"          
## [46] "ASV0210"           "ASV0212"           "ASV0218"          
## [49] "ASV0220"           "ASV0225"           "ASV0262"          
## [52] "ASV0276"           "ASV0289"           "ASV0311"          
## [55] "sum"               "Mesocosm_id"       "Microbe_treatment"
## [58] "Food_type"         "Biomass_day40"
#df.all.larv2 <- df.all.larv1[df.all.larv1$sum!=0,]
#df.all.larv2$sum==0

df.all.larv.long <- melt(df.all.larv1, id.vars=c("sum", "Mesocosm_id", "Microbe_treatment","Food_type", "Biomass_day40"))

df.all.larv.long$rowid <- 1:nrow(df.all.larv.long)

##save
#saveRDS(df.all.water.long,file="df.all.water.long.rds")
#saveRDS(df.all.larv.long,file="df.all.larv.long.rds")

3 Running the statz

Can start from here

setwd("~/nicolagk@hawaii.edu - Google Drive/My Drive/Bamboo_mesos/Bamboo_mesos/04.mixed_models")

df.water.long <- readRDS("df.water.long.rds")
df.larv.long <- readRDS("df.larv.long.rds")

df.all.water.long <- readRDS("df.all.water.long.rds")
df.all.larv.long <- readRDS("df.all.larv.long.rds")

3.1 Water

str(df.water.long)
## 'data.frame':    1860 obs. of  8 variables:
##  $ sum              : num  101694 104508 144484 122042 88536 ...
##  $ Mesocosm_id      : chr  "B_MMO.4" "B_MMO.5" "B_EMO.1" "B_EMO.2" ...
##  $ Microbe_treatment: chr  "MMO" "MMO" "EMO" "EMO" ...
##  $ Food_type        : chr  "Bamboo" "Bamboo" "Bamboo" "Bamboo" ...
##  $ Biomass_day40    : num  2.612 0.468 3.949 0 0.621 ...
##  $ variable         : Factor w/ 62 levels "ASV0002","ASV0003",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ value            : num  26733 27492 5785 9175 5203 ...
##  $ rowid            : int  1 2 3 4 5 6 7 8 9 10 ...
df.water.long$rowid <- as.factor(df.water.long$rowid)
df.water.long$Microbe_treatment <- as.factor(df.water.long$Microbe_treatment)
df.water.long$Food_type <- as.factor(df.water.long$Food_type)
#df.water.long$food.microbes <- as.factor(df.water.long$food.microbes)
df.water.long$Mesocosm_id <- as.factor(df.water.long$Mesocosm_id)
str(df.water.long)
## 'data.frame':    1860 obs. of  8 variables:
##  $ sum              : num  101694 104508 144484 122042 88536 ...
##  $ Mesocosm_id      : Factor w/ 30 levels "B_EMO.1","B_EMO.2",..: 14 15 1 2 3 4 5 6 7 8 ...
##  $ Microbe_treatment: Factor w/ 3 levels "EMO","MEM","MMO": 3 3 1 1 1 1 1 2 2 2 ...
##  $ Food_type        : Factor w/ 2 levels "Bamboo","Larval food": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Biomass_day40    : num  2.612 0.468 3.949 0 0.621 ...
##  $ variable         : Factor w/ 62 levels "ASV0002","ASV0003",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ value            : num  26733 27492 5785 9175 5203 ...
##  $ rowid            : Factor w/ 1860 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
##full model
water.mod.nb1 <- glmmTMB(value~offset(log(sum))+Microbe_treatment*Food_type+(1|variable)+(1|variable:Microbe_treatment)+(1|variable:Food_type)+(1|variable:Microbe_treatment:Food_type),family=nbinom1,dispformula=~variable+Microbe_treatment+Food_type,data=df.water.long)
summary(water.mod.nb1)
##  Family: nbinom1  ( log )
## Formula:          
## value ~ offset(log(sum)) + Microbe_treatment * Food_type + (1 |  
##     variable) + (1 | variable:Microbe_treatment) + (1 | variable:Food_type) +  
##     (1 | variable:Microbe_treatment:Food_type)
## Dispersion:             ~variable + Microbe_treatment + Food_type
## Data: df.water.long
## 
##      AIC      BIC   logLik deviance df.resid 
##  15413.4  15828.1  -7631.7  15263.4     1785 
## 
## Random effects:
## 
## Conditional model:
##  Groups                               Name        Variance  Std.Dev.
##  variable                             (Intercept) 1.377e-06 0.001173
##  variable:Microbe_treatment           (Intercept) 5.699e+00 2.387351
##  variable:Food_type                   (Intercept) 7.351e+00 2.711233
##  variable:Microbe_treatment:Food_type (Intercept) 1.430e-01 0.378103
## Number of obs: 1860, groups:  
## variable, 62; variable:Microbe_treatment, 186; variable:Food_type, 124; variable:Microbe_treatment:Food_type, 372
## 
## Conditional model:
##                                           Estimate Std. Error z value Pr(>|z|)
## (Intercept)                                -6.0382     0.4739 -12.743  < 2e-16
## Microbe_treatmentMEM                       -0.1270     0.4487  -0.283 0.777222
## Microbe_treatmentMMO                       -2.7278     0.6008  -4.540 5.62e-06
## Food_typeLarval food                       -2.0534     0.5393  -3.807 0.000141
## Microbe_treatmentMEM:Food_typeLarval food   0.3195     0.1856   1.722 0.085124
## Microbe_treatmentMMO:Food_typeLarval food  -0.2732     0.3673  -0.744 0.456962
##                                              
## (Intercept)                               ***
## Microbe_treatmentMEM                         
## Microbe_treatmentMMO                      ***
## Food_typeLarval food                      ***
## Microbe_treatmentMEM:Food_typeLarval food .  
## Microbe_treatmentMMO:Food_typeLarval food    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Dispersion model:
##                        Estimate Std. Error z value Pr(>|z|)    
## (Intercept)           6.2564665  0.2975009  21.030  < 2e-16 ***
## variableASV0003      -1.5088503  0.5166781  -2.920 0.003497 ** 
## variableASV0004       1.8987349  0.4512348   4.208 2.58e-05 ***
## variableASV0005       1.8317243  0.4502982   4.068 4.75e-05 ***
## variableASV0007       0.5783152  0.4281055   1.351 0.176737    
## variableASV0008       0.2901073  0.4510498   0.643 0.520106    
## variableASV0009      -0.1788562  0.4028234  -0.444 0.657038    
## variableASV0010       0.7611656  0.4562629   1.668 0.095264 .  
## variableASV0011       1.3035506  0.4568341   2.853 0.004325 ** 
## variableASV0012       1.2093853  0.4284904   2.822 0.004766 ** 
## variableASV0015       0.8622593  0.4869641   1.771 0.076613 .  
## variableASV0016      -0.2424918  0.4650040  -0.521 0.602030    
## variableASV0017      -0.5906261  0.5638137  -1.048 0.294843    
## variableASV0020       1.3479492  0.4529875   2.976 0.002923 ** 
## variableASV0021       0.1122259  0.7669599   0.146 0.883664    
## variableASV0022       0.5214020  0.4492747   1.161 0.245828    
## variableASV0023       0.6708045  0.6555633   1.023 0.306190    
## variableASV0024       1.2393792  0.4658077   2.661 0.007798 ** 
## variableASV0025      -0.3619200  0.6020512  -0.601 0.547743    
## variableASV0027      -1.0459947  0.5842380  -1.790 0.073397 .  
## variableASV0036       3.6861632  0.7364961   5.005 5.59e-07 ***
## variableASV0037       2.6200524  0.9447735   2.773 0.005551 ** 
## variableASV0039       0.0883535  0.5842289   0.151 0.879793    
## variableASV0040      -0.5393587  0.6602007  -0.817 0.413950    
## variableASV0042       1.8028771  0.5040531   3.577 0.000348 ***
## variableASV0044       0.8508108  0.5877397   1.448 0.147730    
## variableASV0051      -1.1642234  0.4504576  -2.585 0.009751 ** 
## variableASV0054      -0.3722326  0.6419720  -0.580 0.562031    
## variableASV0057      -0.6032034  0.5873371  -1.027 0.304414    
## variableASV0059      -0.2193775  0.5489210  -0.400 0.689413    
## variableASV0063      -0.4682974  0.5026902  -0.932 0.351552    
## variableASV0065      -0.4523079  0.6357944  -0.711 0.476833    
## variableASV0066       0.9424386  0.7085647   1.330 0.183496    
## variableASV0069      -1.9006165  0.6244191  -3.044 0.002336 ** 
## variableASV0072      -1.7610774  0.6163257  -2.857 0.004272 ** 
## variableASV0073       1.6480567  0.7696031   2.141 0.032239 *  
## variableASV0079      -1.3420093  0.7566616  -1.774 0.076131 .  
## variableASV0080       1.0856345  0.6405383   1.695 0.090099 .  
## variableASV0082       0.4498012  0.5686437   0.791 0.428940    
## variableASV0084      -0.2608202  0.5945235  -0.439 0.660876    
## variableASV0085      -1.1383417  0.4505072  -2.527 0.011511 *  
## variableASV0089       0.0001207  0.6510603   0.000 0.999852    
## variableASV0091       0.2946653  0.7864269   0.375 0.707892    
## variableASV0093      -1.8915201  0.6102519  -3.100 0.001938 ** 
## variableASV0094      -0.5932495  0.6157738  -0.963 0.335336    
## variableASV0095       0.4885778  0.6624257   0.738 0.460783    
## variableASV0101       1.2566686  0.8443749   1.488 0.136676    
## variableASV0102      -0.3685071  0.6095838  -0.605 0.545496    
## variableASV0105       0.8932714  0.9708123   0.920 0.357506    
## variableASV0111      -1.4285033  0.5539093  -2.579 0.009910 ** 
## variableASV0112       0.7922081  0.7841482   1.010 0.312362    
## variableASV0115       0.0603164  0.6428722   0.094 0.925250    
## variableASV0116      -0.1551065  0.5273773  -0.294 0.768674    
## variableASV0119      -0.2673950  0.6575854  -0.407 0.684279    
## variableASV0123      -0.6364850  0.6579777  -0.967 0.333376    
## variableASV0128       0.7760974  0.9009622   0.861 0.389013    
## variableASV0139       1.1327518  0.8767022   1.292 0.196336    
## variableASV0145      -2.6420271  0.8629926  -3.061 0.002203 ** 
## variableASV0146       0.0246931  0.7256392   0.034 0.972854    
## variableASV0155      -1.9592978  0.6293727  -3.113 0.001851 ** 
## variableASV0167      -0.5380176  0.7807569  -0.689 0.490762    
## variableASV0180      -3.2500265  0.6435964  -5.050 4.42e-07 ***
## Microbe_treatmentMEM  0.0107150  0.1390419   0.077 0.938573    
## Microbe_treatmentMMO  2.0624551  0.2556546   8.067 7.18e-16 ***
## Food_typeLarval food -0.1704443  0.1641874  -1.038 0.299220    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
water.mod.nb2 <- glmmTMB(value~offset(log(sum))+Microbe_treatment*Food_type+(1|variable)+(1|variable:Microbe_treatment)+(1|variable:Food_type)+(1|variable:Microbe_treatment:Food_type),family=nbinom2,dispformula=~variable+Microbe_treatment+Food_type,data=df.water.long)
summary(water.mod.nb2)
##  Family: nbinom2  ( log )
## Formula:          
## value ~ offset(log(sum)) + Microbe_treatment * Food_type + (1 |  
##     variable) + (1 | variable:Microbe_treatment) + (1 | variable:Food_type) +  
##     (1 | variable:Microbe_treatment:Food_type)
## Dispersion:             ~variable + Microbe_treatment + Food_type
## Data: df.water.long
## 
##      AIC      BIC   logLik deviance df.resid 
##  15930.5  16345.2  -7890.3  15780.5     1785 
## 
## Random effects:
## 
## Conditional model:
##  Groups                               Name        Variance  Std.Dev. 
##  variable                             (Intercept) 9.984e-07 0.0009992
##  variable:Microbe_treatment           (Intercept) 2.456e+01 4.9561702
##  variable:Food_type                   (Intercept) 2.132e+01 4.6173294
##  variable:Microbe_treatment:Food_type (Intercept) 1.057e-01 0.3250761
## Number of obs: 1860, groups:  
## variable, 62; variable:Microbe_treatment, 186; variable:Food_type, 124; variable:Microbe_treatment:Food_type, 372
## 
## Conditional model:
##                                           Estimate Std. Error z value Pr(>|z|)
## (Intercept)                               -6.35573    0.87820  -7.237 4.58e-13
## Microbe_treatmentMEM                      -0.05179    0.91574  -0.057   0.9549
## Microbe_treatmentMMO                      -8.93860    1.69505  -5.273 1.34e-07
## Food_typeLarval food                      -3.67418    0.90354  -4.066 4.77e-05
## Microbe_treatmentMEM:Food_typeLarval food  0.30694    0.20736   1.480   0.1388
## Microbe_treatmentMMO:Food_typeLarval food -0.74827    0.43617  -1.716   0.0862
##                                              
## (Intercept)                               ***
## Microbe_treatmentMEM                         
## Microbe_treatmentMMO                      ***
## Food_typeLarval food                      ***
## Microbe_treatmentMEM:Food_typeLarval food    
## Microbe_treatmentMMO:Food_typeLarval food .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Dispersion model:
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)           2.44408    0.28413   8.602  < 2e-16 ***
## variableASV0003       1.01324    0.46818   2.164 0.030447 *  
## variableASV0004      -0.83267    0.43889  -1.897 0.057801 .  
## variableASV0005      -0.59799    0.43798  -1.365 0.172146    
## variableASV0007      -1.29878    0.39888  -3.256 0.001130 ** 
## variableASV0008      -2.43096    0.43906  -5.537 3.08e-08 ***
## variableASV0009      -0.17048    0.38863  -0.439 0.660899    
## variableASV0010      -0.53041    0.44308  -1.197 0.231267    
## variableASV0011      -1.55520    0.42533  -3.656 0.000256 ***
## variableASV0012      -2.25799    0.40576  -5.565 2.62e-08 ***
## variableASV0015      -2.91623    0.41250  -7.070 1.55e-12 ***
## variableASV0016       0.63873    0.45136   1.415 0.157030    
## variableASV0017      -2.22143    0.52617  -4.222 2.42e-05 ***
## variableASV0020      -1.77903    0.41558  -4.281 1.86e-05 ***
## variableASV0021      -3.33457    0.55482  -6.010 1.85e-09 ***
## variableASV0022      -1.72471    0.41958  -4.111 3.95e-05 ***
## variableASV0023      -2.88667    0.54278  -5.318 1.05e-07 ***
## variableASV0024      -2.23507    0.42275  -5.287 1.24e-07 ***
## variableASV0025      -0.65671    0.55776  -1.177 0.239035    
## variableASV0027      -2.19496    0.55443  -3.959 7.53e-05 ***
## variableASV0036      -4.71917    0.44817 -10.530  < 2e-16 ***
## variableASV0037      -5.78499    0.52382 -11.044  < 2e-16 ***
## variableASV0039      -3.91863    0.49526  -7.912 2.53e-15 ***
## variableASV0040      -2.58160    0.55868  -4.621 3.82e-06 ***
## variableASV0042      -3.02878    0.41361  -7.323 2.43e-13 ***
## variableASV0044      -1.55150    0.53856  -2.881 0.003966 ** 
## variableASV0051      -2.43832    0.43082  -5.660 1.52e-08 ***
## variableASV0054      -3.98059    0.45992  -8.655  < 2e-16 ***
## variableASV0057      -3.48125    0.48176  -7.226 4.97e-13 ***
## variableASV0059      -3.74888    0.45429  -8.252  < 2e-16 ***
## variableASV0063      -2.99035    0.44526  -6.716 1.87e-11 ***
## variableASV0065      -3.32762    0.51256  -6.492 8.46e-11 ***
## variableASV0066      -4.30290    0.48849  -8.809  < 2e-16 ***
## variableASV0069       0.27072    0.57313   0.472 0.636674    
## variableASV0072      -0.13073    0.56912  -0.230 0.818321    
## variableASV0073      -4.73197    0.49269  -9.604  < 2e-16 ***
## variableASV0079      -2.66730    0.57854  -4.610 4.02e-06 ***
## variableASV0080      -4.13395    0.45285  -9.129  < 2e-16 ***
## variableASV0082      -3.69105    0.43927  -8.403  < 2e-16 ***
## variableASV0084      -3.01801    0.47206  -6.393 1.62e-10 ***
## variableASV0085      -1.19092    0.42990  -2.770 0.005602 ** 
## variableASV0089      -2.62620    0.55567  -4.726 2.29e-06 ***
## variableASV0091      -3.87883    0.57103  -6.793 1.10e-11 ***
## variableASV0093      -0.61818    0.56392  -1.096 0.272979    
## variableASV0094      -0.64417    0.56710  -1.136 0.256000    
## variableASV0095      -2.91679    0.53455  -5.457 4.85e-08 ***
## variableASV0101      -4.77505    0.49073  -9.731  < 2e-16 ***
## variableASV0102      -2.95578    0.48183  -6.135 8.54e-10 ***
## variableASV0105      -4.40262    0.61297  -7.182 6.85e-13 ***
## variableASV0111      -3.66298    0.46197  -7.929 2.21e-15 ***
## variableASV0112      -3.52209    0.54224  -6.495 8.28e-11 ***
## variableASV0115      -4.23469    0.49012  -8.640  < 2e-16 ***
## variableASV0116      -3.54838    0.42923  -8.267  < 2e-16 ***
## variableASV0119      -2.56979    0.56030  -4.586 4.51e-06 ***
## variableASV0123      -2.68720    0.55503  -4.842 1.29e-06 ***
## variableASV0128      -3.90753    0.56794  -6.880 5.98e-12 ***
## variableASV0139      -3.97281    0.56553  -7.025 2.14e-12 ***
## variableASV0145      -2.14679    0.61828  -3.472 0.000516 ***
## variableASV0146      -4.27875    0.48611  -8.802  < 2e-16 ***
## variableASV0155      -0.09976    0.57780  -0.173 0.862922    
## variableASV0167      -3.60614    0.58520  -6.162 7.17e-10 ***
## variableASV0180       0.16370    0.59327   0.276 0.782605    
## Microbe_treatmentMEM -0.18359    0.12110  -1.516 0.129514    
## Microbe_treatmentMMO -1.38803    0.25430  -5.458 4.81e-08 ***
## Food_typeLarval food -0.20719    0.14201  -1.459 0.144567    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
AICtab(water.mod.nb1,water.mod.nb2)
##               dAIC  df
## water.mod.nb1   0.0 75
## water.mod.nb2 517.1 75
water.mod.nb1.resid <- simulateResiduals(fittedModel = water.mod.nb1, plot = T)

water.mod.nb2.resid <- simulateResiduals(fittedModel = water.mod.nb2, plot = T)

##dispersal formula checks
water.mod.nb1.disnovar <- glmmTMB(value~offset(log(sum))+Microbe_treatment*Food_type+(1|variable)+(1|variable:Microbe_treatment)+(1|variable:Food_type)+(1|variable:Microbe_treatment:Food_type),family=nbinom1,dispformula=~Microbe_treatment+Food_type,data=df.water.long)

anova(water.mod.nb1.disnovar,water.mod.nb1) #sig***
## Data: df.water.long
## Models:
## water.mod.nb1.disnovar: value ~ offset(log(sum)) + Microbe_treatment * Food_type + (1 | , zi=~0, disp=~Microbe_treatment + Food_type
## water.mod.nb1.disnovar:     variable) + (1 | variable:Microbe_treatment) + (1 | variable:Food_type) + , zi=~0, disp=~variable + Microbe_treatment + Food_type
## water.mod.nb1.disnovar:     (1 | variable:Microbe_treatment:Food_type), zi=~0, disp=~Microbe_treatment + Food_type
## water.mod.nb1: value ~ offset(log(sum)) + Microbe_treatment * Food_type + (1 | , zi=~0, disp=~variable + Microbe_treatment + Food_type
## water.mod.nb1:     variable) + (1 | variable:Microbe_treatment) + (1 | variable:Food_type) + , zi=~0, disp=~Microbe_treatment + Food_type
## water.mod.nb1:     (1 | variable:Microbe_treatment:Food_type), zi=~0, disp=~variable + Microbe_treatment + Food_type
##                        Df   AIC   BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## water.mod.nb1.disnovar 14 15642 15720 -7807.1    15614                         
## water.mod.nb1          75 15413 15828 -7631.7    15263 350.86     61  < 2.2e-16
##                           
## water.mod.nb1.disnovar    
## water.mod.nb1          ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##removing 3 way interaction
water.mod.nb1.nofm <- update(water.mod.nb1,.~.-(1|variable:Microbe_treatment:Food_type))
anova(water.mod.nb1,water.mod.nb1.nofm)
## Data: df.water.long
## Models:
## water.mod.nb1.nofm: value ~ Microbe_treatment + Food_type + (1 | variable) + (1 | , zi=~0, disp=~variable + Microbe_treatment + Food_type
## water.mod.nb1.nofm:     variable:Microbe_treatment) + (1 | variable:Food_type) + , zi=~0, disp=~variable + Microbe_treatment + Food_type
## water.mod.nb1.nofm:     Microbe_treatment:Food_type + offset(log(sum)), zi=~0, disp=~variable + Microbe_treatment + Food_type
## water.mod.nb1: value ~ offset(log(sum)) + Microbe_treatment * Food_type + (1 | , zi=~0, disp=~variable + Microbe_treatment + Food_type
## water.mod.nb1:     variable) + (1 | variable:Microbe_treatment) + (1 | variable:Food_type) + , zi=~0, disp=~variable + Microbe_treatment + Food_type
## water.mod.nb1:     (1 | variable:Microbe_treatment:Food_type), zi=~0, disp=~variable + Microbe_treatment + Food_type
##                    Df   AIC   BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
## water.mod.nb1.nofm 74 15425 15834 -7638.4    15277                             
## water.mod.nb1      75 15413 15828 -7631.7    15263 13.473      1   0.000242 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##super sig 0.000242

water.mod.nb2.nofm <- update(water.mod.nb2,.~.-(1|variable:Microbe_treatment:Food_type))
anova(water.mod.nb2,water.mod.nb2.nofm)
## Data: df.water.long
## Models:
## water.mod.nb2.nofm: value ~ Microbe_treatment + Food_type + (1 | variable) + (1 | , zi=~0, disp=~variable + Microbe_treatment + Food_type
## water.mod.nb2.nofm:     variable:Microbe_treatment) + (1 | variable:Food_type) + , zi=~0, disp=~variable + Microbe_treatment + Food_type
## water.mod.nb2.nofm:     Microbe_treatment:Food_type + offset(log(sum)), zi=~0, disp=~variable + Microbe_treatment + Food_type
## water.mod.nb2: value ~ offset(log(sum)) + Microbe_treatment * Food_type + (1 | , zi=~0, disp=~variable + Microbe_treatment + Food_type
## water.mod.nb2:     variable) + (1 | variable:Microbe_treatment) + (1 | variable:Food_type) + , zi=~0, disp=~variable + Microbe_treatment + Food_type
## water.mod.nb2:     (1 | variable:Microbe_treatment:Food_type), zi=~0, disp=~variable + Microbe_treatment + Food_type
##                    Df   AIC   BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)   
## water.mod.nb2.nofm 74 15937 16346 -7894.4    15789                            
## water.mod.nb2      75 15930 16345 -7890.3    15780 8.3132      1   0.003936 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##also super sig

###save conditional modes
water.mod.ranef <- as.data.frame(ranef(water.mod.nb1))
#write.csv(water.mod.ranef, "water.mod.ranef.csv")
#saveRDS(water.mod.nb1,file="water.mod.nb1.rds")

3.1.1 Nbinom1 vs. nbinom2 things

Following this wonderful resource, on page 16

Ben Bolker, Mollie Brooks, Beth Gardner, Cleridy Lennert, Mihoko Minami, October 23, 2012, Owls example: a zero-inflated, generalized linear mixed model for count data. https://groups.nceas.ucsb.edu/non-linear-modeling/projects/owls/WRITEUP/owls.pdf

mvtab <- ddply(df.water.long,
               .(variable:Microbe_treatment:Food_type),
               summarise,
               callmean=mean(value),
               callvar=var(value))

q1 <- qplot(callmean,callvar,data=mvtab)
## Warning: `qplot()` was deprecated in ggplot2 3.4.0.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
print(q1+
        ## linear (quasi-Poisson/NB1) fit
        geom_smooth(method="lm",formula=y~x-1,colour="pink")+
        ## smooth (loess)
        geom_smooth(colour="red")+
        ## semi-quadratic (NB2/LNP)
        geom_smooth(method="lm",formula=y~I(x^2)+offset(x)-1,colour="purple")+
        ## Poisson (v=m)
        geom_abline(a=0,b=1,lty=2))
## Warning in geom_abline(a = 0, b = 1, lty = 2): Ignoring unknown parameters: `a`
## and `b`
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

3.2 Larvae

str(df.larv.long)
## 'data.frame':    1320 obs. of  8 variables:
##  $ sum              : num  63543 18211 13437 48065 39128 ...
##  $ Mesocosm_id      : chr  "B_EMO.3" "B_EMO.4" "B_EMO.5" "B_MEM.1" ...
##  $ Microbe_treatment: chr  "EMO" "EMO" "EMO" "MEM" ...
##  $ Food_type        : chr  "Bamboo" "Bamboo" "Bamboo" "Bamboo" ...
##  $ Biomass_day40    : num  0.621 0.308 4.085 0.227 4.007 ...
##  $ variable         : Factor w/ 44 levels "ASV0002","ASV0003",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ value            : num  0 20 0 0 0 42 0 0 257 0 ...
##  $ rowid            : int  1 2 3 4 5 6 7 8 9 10 ...
df.larv.long$rowid <- as.factor(df.larv.long$rowid)
df.larv.long$Microbe_treatment <- as.factor(df.larv.long$Microbe_treatment)
df.larv.long$Food_type <- as.factor(df.larv.long$Food_type)
#df.larv.long$food.microbes <- as.factor(df.larv.long$food.microbes)
df.larv.long$Mesocosm_id <- as.factor(df.larv.long$Mesocosm_id)
str(df.larv.long)
## 'data.frame':    1320 obs. of  8 variables:
##  $ sum              : num  63543 18211 13437 48065 39128 ...
##  $ Mesocosm_id      : Factor w/ 30 levels "B_EMO.1","B_EMO.2",..: 3 4 5 6 7 8 9 10 11 12 ...
##  $ Microbe_treatment: Factor w/ 3 levels "EMO","MEM","MMO": 1 1 1 2 2 2 2 2 3 3 ...
##  $ Food_type        : Factor w/ 2 levels "Bamboo","Larval food": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Biomass_day40    : num  0.621 0.308 4.085 0.227 4.007 ...
##  $ variable         : Factor w/ 44 levels "ASV0002","ASV0003",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ value            : num  0 20 0 0 0 42 0 0 257 0 ...
##  $ rowid            : Factor w/ 1320 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
##binomial
larv.mod.bin <- glmmTMB(cbind(value, sum-value)~Microbe_treatment*Food_type+(1|variable)+(1|variable:Microbe_treatment)+(1|variable:Food_type)+(1|variable:Microbe_treatment:Food_type)+(1|rowid),family=binomial,data=df.larv.long)
summary(larv.mod.bin)
##  Family: binomial  ( logit )
## Formula:          
## cbind(value, sum - value) ~ Microbe_treatment * Food_type + (1 |  
##     variable) + (1 | variable:Microbe_treatment) + (1 | variable:Food_type) +  
##     (1 | variable:Microbe_treatment:Food_type) + (1 | rowid)
## Data: df.larv.long
## 
##      AIC      BIC   logLik deviance df.resid 
##   9953.6  10010.7  -4965.8   9931.6     1309 
## 
## Random effects:
## 
## Conditional model:
##  Groups                               Name        Variance  Std.Dev. 
##  variable                             (Intercept) 6.837e-07 0.0008269
##  variable:Microbe_treatment           (Intercept) 2.057e+01 4.5355850
##  variable:Food_type                   (Intercept) 3.687e+01 6.0717861
##  variable:Microbe_treatment:Food_type (Intercept) 1.325e+01 3.6401214
##  rowid                                (Intercept) 1.003e+01 3.1666214
## Number of obs: 1320, groups:  
## variable, 44; variable:Microbe_treatment, 132; variable:Food_type, 88; variable:Microbe_treatment:Food_type, 264; rowid, 1320
## 
## Conditional model:
##                                             Estimate Std. Error z value
## (Intercept)                               -14.040307   1.440472  -9.747
## Microbe_treatmentMEM                        3.377608   1.410992   2.394
## Microbe_treatmentMMO                       -4.446516   1.613997  -2.755
## Food_typeLarval food                        0.003874   1.729333   0.002
## Microbe_treatmentMEM:Food_typeLarval food  -3.474420   1.498526  -2.319
## Microbe_treatmentMMO:Food_typeLarval food  -2.067792   1.602383  -1.290
##                                           Pr(>|z|)    
## (Intercept)                                < 2e-16 ***
## Microbe_treatmentMEM                       0.01668 *  
## Microbe_treatmentMMO                       0.00587 ** 
## Food_typeLarval food                       0.99821    
## Microbe_treatmentMEM:Food_typeLarval food  0.02042 *  
## Microbe_treatmentMMO:Food_typeLarval food  0.19689    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##full model
larv.mod.nb1 <- glmmTMB(value~offset(log(sum))+Microbe_treatment*Food_type+(1|variable)+(1|variable:Microbe_treatment)+(1|variable:Food_type)+(1|variable:Microbe_treatment:Food_type),disp=~variable+Food_type+Microbe_treatment,family=nbinom1,data=df.larv.long)
summary(larv.mod.nb1)
##  Family: nbinom1  ( log )
## Formula:          
## value ~ offset(log(sum)) + Microbe_treatment * Food_type + (1 |  
##     variable) + (1 | variable:Microbe_treatment) + (1 | variable:Food_type) +  
##     (1 | variable:Microbe_treatment:Food_type)
## Dispersion:             ~variable + Food_type + Microbe_treatment
## Data: df.larv.long
## 
##      AIC      BIC   logLik deviance df.resid 
##   9374.9   9670.5  -4630.5   9260.9     1263 
## 
## Random effects:
## 
## Conditional model:
##  Groups                               Name        Variance  Std.Dev. 
##  variable                             (Intercept) 1.856e-07 0.0004308
##  variable:Microbe_treatment           (Intercept) 3.439e+00 1.8544060
##  variable:Food_type                   (Intercept) 7.217e+00 2.6865139
##  variable:Microbe_treatment:Food_type (Intercept) 6.505e-01 0.8065311
## Number of obs: 1320, groups:  
## variable, 44; variable:Microbe_treatment, 132; variable:Food_type, 88; variable:Microbe_treatment:Food_type, 264
## 
## Conditional model:
##                                           Estimate Std. Error z value Pr(>|z|)
## (Intercept)                               -6.95183    0.58047 -11.976  < 2e-16
## Microbe_treatmentMEM                       0.29763    0.50828   0.586  0.55817
## Microbe_treatmentMMO                      -1.62053    0.59085  -2.743  0.00609
## Food_typeLarval food                       0.17975    0.70299   0.256  0.79819
## Microbe_treatmentMEM:Food_typeLarval food -0.53532    0.42755  -1.252  0.21055
## Microbe_treatmentMMO:Food_typeLarval food -0.06406    0.52157  -0.123  0.90225
##                                              
## (Intercept)                               ***
## Microbe_treatmentMEM                         
## Microbe_treatmentMMO                      ** 
## Food_typeLarval food                         
## Microbe_treatmentMEM:Food_typeLarval food    
## Microbe_treatmentMMO:Food_typeLarval food    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Dispersion model:
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)           5.56954    0.58181   9.573  < 2e-16 ***
## variableASV0003       2.42329    0.62861   3.855 0.000116 ***
## variableASV0005       1.09329    0.85231   1.283 0.199582    
## variableASV0006       2.70004    0.65837   4.101 4.11e-05 ***
## variableASV0007       1.55254    0.66242   2.344 0.019092 *  
## variableASV0009       1.99002    0.69265   2.873 0.004065 ** 
## variableASV0010       0.91182    0.73202   1.246 0.212905    
## variableASV0011      -0.13271    0.75148  -0.177 0.859820    
## variableASV0012       1.80874    1.01104   1.789 0.073615 .  
## variableASV0013       2.91827    0.76372   3.821 0.000133 ***
## variableASV0016       1.53602    0.85862   1.789 0.073625 .  
## variableASV0019       2.82165    0.71026   3.973 7.11e-05 ***
## variableASV0024       1.29458    0.70167   1.845 0.065037 .  
## variableASV0027       1.97844    0.72858   2.715 0.006618 ** 
## variableASV0032       2.00643    0.79047   2.538 0.011140 *  
## variableASV0033       2.64198    0.71766   3.681 0.000232 ***
## variableASV0043       4.43006    0.87997   5.034 4.79e-07 ***
## variableASV0053       1.16523    0.72785   1.601 0.109397    
## variableASV0054       1.68061    0.89688   1.874 0.060952 .  
## variableASV0055       2.38325    0.78221   3.047 0.002313 ** 
## variableASV0058       2.93112    0.84328   3.476 0.000509 ***
## variableASV0059       0.60771    0.70774   0.859 0.390530    
## variableASV0060       3.10043    0.83244   3.725 0.000196 ***
## variableASV0061       1.71121    0.81456   2.101 0.035660 *  
## variableASV0066      -0.47107    0.73958  -0.637 0.524162    
## variableASV0070       1.37715    0.92483   1.489 0.136467    
## variableASV0079       0.76562    0.68258   1.122 0.262005    
## variableASV0083       0.62041    0.68273   0.909 0.363494    
## variableASV0085       1.21442    1.03035   1.179 0.238536    
## variableASV0087       1.81063    0.76428   2.369 0.017833 *  
## variableASV0099       1.67867    0.82274   2.040 0.041317 *  
## variableASV0100       1.00298    0.74743   1.342 0.179630    
## variableASV0102      -0.43747    0.85386  -0.512 0.608412    
## variableASV0125       2.27182    0.87847   2.586 0.009707 ** 
## variableASV0130       0.98610    0.71145   1.386 0.165730    
## variableASV0134       2.22981    1.00464   2.220 0.026452 *  
## variableASV0138       0.41231    0.88483   0.466 0.641228    
## variableASV0149      -1.11486    0.70938  -1.572 0.116043    
## variableASV0150       0.06528    0.81195   0.080 0.935924    
## variableASV0164      -0.18252    0.73054  -0.250 0.802708    
## variableASV0218       0.10549    0.94527   0.112 0.911145    
## variableASV0262       0.74410    1.03491   0.719 0.472138    
## variableASV0289      -0.15667    0.86232  -0.182 0.855832    
## variableASV0311       0.23946    0.88198   0.272 0.786003    
## Food_typeLarval food  0.93604    0.28051   3.337 0.000847 ***
## Microbe_treatmentMEM -0.77455    0.23725  -3.265 0.001096 ** 
## Microbe_treatmentMMO  0.05033    0.31572   0.159 0.873344    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#saveRDS(larv.mod.nb1,file="larv.mod.nb1.rds")

larv.mod.nb2 <- glmmTMB(value~offset(log(sum))+Microbe_treatment*Food_type+(1|variable)+(1|variable:Microbe_treatment)+(1|variable:Food_type)+(1|variable:Microbe_treatment:Food_type),disp=~variable+Food_type+Microbe_treatment,family=nbinom2,data=df.larv.long)
summary(larv.mod.nb2)
##  Family: nbinom2  ( log )
## Formula:          
## value ~ offset(log(sum)) + Microbe_treatment * Food_type + (1 |  
##     variable) + (1 | variable:Microbe_treatment) + (1 | variable:Food_type) +  
##     (1 | variable:Microbe_treatment:Food_type)
## Dispersion:             ~variable + Food_type + Microbe_treatment
## Data: df.larv.long
## 
##      AIC      BIC   logLik deviance df.resid 
##   9752.6  10048.2  -4819.3   9638.6     1263 
## 
## Random effects:
## 
## Conditional model:
##  Groups                               Name        Variance  Std.Dev.
##  variable                             (Intercept) 7.414e-07 0.000861
##  variable:Microbe_treatment           (Intercept) 1.865e+01 4.318578
##  variable:Food_type                   (Intercept) 2.931e+01 5.413738
##  variable:Microbe_treatment:Food_type (Intercept) 3.068e+00 1.751673
## Number of obs: 1320, groups:  
## variable, 44; variable:Microbe_treatment, 132; variable:Food_type, 88; variable:Microbe_treatment:Food_type, 264
## 
## Conditional model:
##                                           Estimate Std. Error z value Pr(>|z|)
## (Intercept)                               -9.53380    1.18259  -8.062 7.52e-16
## Microbe_treatmentMEM                       2.60767    1.11997   2.328 0.019894
## Microbe_treatmentMMO                      -2.67102    1.25956  -2.121 0.033955
## Food_typeLarval food                       0.01895    1.37397   0.014 0.988998
## Microbe_treatmentMEM:Food_typeLarval food -2.99036    0.89874  -3.327 0.000877
## Microbe_treatmentMMO:Food_typeLarval food -2.69280    1.10875  -2.429 0.015154
##                                              
## (Intercept)                               ***
## Microbe_treatmentMEM                      *  
## Microbe_treatmentMMO                      *  
## Food_typeLarval food                         
## Microbe_treatmentMEM:Food_typeLarval food ***
## Microbe_treatmentMMO:Food_typeLarval food *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Dispersion model:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)          -2.779935   0.372728  -7.458 8.76e-14 ***
## variableASV0003       1.636689   0.442891   3.695 0.000219 ***
## variableASV0005       0.518260   0.591254   0.877 0.380735    
## variableASV0006       1.595294   0.496609   3.212 0.001316 ** 
## variableASV0007       1.145117   0.485977   2.356 0.018457 *  
## variableASV0009       0.830116   0.478662   1.734 0.082875 .  
## variableASV0010       0.653609   0.485088   1.347 0.177851    
## variableASV0011       2.628128   0.585580   4.488 7.19e-06 ***
## variableASV0012      -0.861885   0.648264  -1.330 0.183674    
## variableASV0013       4.181151   0.598683   6.984 2.87e-12 ***
## variableASV0016       0.487424   0.590490   0.825 0.409112    
## variableASV0019       1.813278   0.519091   3.493 0.000477 ***
## variableASV0024       0.947252   0.494014   1.917 0.055179 .  
## variableASV0027       1.446153   0.526480   2.747 0.006017 ** 
## variableASV0032       3.822804   0.592276   6.454 1.09e-10 ***
## variableASV0033       0.702789   0.462448   1.520 0.128583    
## variableASV0043       0.743267   0.610073   1.218 0.223100    
## variableASV0053       0.945001   0.503333   1.877 0.060451 .  
## variableASV0054      -0.804134   0.527172  -1.525 0.127166    
## variableASV0055       2.566564   0.558195   4.598 4.27e-06 ***
## variableASV0058       1.378698   0.593636   2.322 0.020208 *  
## variableASV0059       0.975630   0.512528   1.904 0.056967 .  
## variableASV0060       0.513368   0.536041   0.958 0.338213    
## variableASV0061       0.474951   0.499538   0.951 0.341716    
## variableASV0066       1.136746   0.593298   1.916 0.055368 .  
## variableASV0070       1.898357   0.628794   3.019 0.002536 ** 
## variableASV0079       1.510157   0.478858   3.154 0.001612 ** 
## variableASV0083       1.551743   0.480928   3.227 0.001253 ** 
## variableASV0085      -0.545332   0.610078  -0.894 0.371390    
## variableASV0087       2.811920   0.564989   4.977 6.46e-07 ***
## variableASV0099      -0.396072   0.495977  -0.799 0.424541    
## variableASV0100       0.979286   0.547645   1.788 0.073747 .  
## variableASV0102      -0.494742   0.528771  -0.936 0.349455    
## variableASV0125       0.788877   0.547143   1.442 0.149356    
## variableASV0130       0.169561   0.454711   0.373 0.709224    
## variableASV0134      -0.822009   0.552742  -1.487 0.136976    
## variableASV0138       1.647222   0.599422   2.748 0.005996 ** 
## variableASV0149       3.139577   0.557746   5.629 1.81e-08 ***
## variableASV0150       1.208037   0.582530   2.074 0.038100 *  
## variableASV0164       0.800755   0.511426   1.566 0.117412    
## variableASV0218       0.323392   0.630189   0.513 0.607835    
## variableASV0262       0.118794   0.672205   0.177 0.859726    
## variableASV0289       1.829776   0.615928   2.971 0.002971 ** 
## variableASV0311       1.711515   0.605002   2.829 0.004670 ** 
## Food_typeLarval food  0.884808   0.167399   5.286 1.25e-07 ***
## Microbe_treatmentMEM  0.251416   0.168605   1.491 0.135922    
## Microbe_treatmentMMO -0.008739   0.223269  -0.039 0.968779    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
AICtab(larv.mod.nb1,larv.mod.nb2,larv.mod.bin)
##              dAIC  df
## larv.mod.nb1   0.0 57
## larv.mod.nb2 377.7 57
## larv.mod.bin 578.7 11
larv.mod.bin.resid <- simulateResiduals(fittedModel = larv.mod.bin, plot = T)
## DHARMa:testOutliers with type = binomial may have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details

larv.mod.nb1.resid <- simulateResiduals(fittedModel = larv.mod.nb1, plot = T)
## DHARMa:testOutliers with type = binomial may have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details

larv.mod.nb2.resid <- simulateResiduals(fittedModel = larv.mod.nb2, plot = T)
## DHARMa:testOutliers with type = binomial may have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details

##dispersal formula checks
larv.mod.nb1.disnovar <- glmmTMB(value~offset(log(sum))+Microbe_treatment*Food_type+(1|variable)+(1|variable:Microbe_treatment)+(1|variable:Food_type)+(1|variable:Microbe_treatment:Food_type),family=nbinom1,dispformula=~Microbe_treatment+Food_type,data=df.larv.long)

anova(larv.mod.nb1.disnovar,larv.mod.nb1) #sig***
## Data: df.larv.long
## Models:
## larv.mod.nb1.disnovar: value ~ offset(log(sum)) + Microbe_treatment * Food_type + (1 | , zi=~0, disp=~Microbe_treatment + Food_type
## larv.mod.nb1.disnovar:     variable) + (1 | variable:Microbe_treatment) + (1 | variable:Food_type) + , zi=~0, disp=~variable + Food_type + Microbe_treatment
## larv.mod.nb1.disnovar:     (1 | variable:Microbe_treatment:Food_type), zi=~0, disp=~Microbe_treatment + Food_type
## larv.mod.nb1: value ~ offset(log(sum)) + Microbe_treatment * Food_type + (1 | , zi=~0, disp=~variable + Food_type + Microbe_treatment
## larv.mod.nb1:     variable) + (1 | variable:Microbe_treatment) + (1 | variable:Food_type) + , zi=~0, disp=~Microbe_treatment + Food_type
## larv.mod.nb1:     (1 | variable:Microbe_treatment:Food_type), zi=~0, disp=~variable + Food_type + Microbe_treatment
##                       Df    AIC    BIC  logLik deviance  Chisq Chi Df
## larv.mod.nb1.disnovar 14 9456.8 9529.4 -4714.4   9428.8              
## larv.mod.nb1          57 9374.9 9670.5 -4630.5   9260.9 167.92     43
##                       Pr(>Chisq)    
## larv.mod.nb1.disnovar               
## larv.mod.nb1           < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##removing 3 way interaction
larv.mod.nb1.nofm <- update(larv.mod.nb1,.~.-(1|variable:Microbe_treatment:Food_type))
summary(larv.mod.nb1.nofm)
##  Family: nbinom1  ( log )
## Formula:          
## value ~ Microbe_treatment + Food_type + (1 | variable) + (1 |  
##     variable:Microbe_treatment) + (1 | variable:Food_type) +  
##     Microbe_treatment:Food_type + offset(log(sum))
## Dispersion:             ~variable + Food_type + Microbe_treatment
## Data: df.larv.long
## 
##      AIC      BIC   logLik deviance df.resid 
##   9380.4   9670.8  -4634.2   9268.4     1264 
## 
## Random effects:
## 
## Conditional model:
##  Groups                     Name        Variance  Std.Dev. 
##  variable                   (Intercept) 1.833e-07 0.0004282
##  variable:Microbe_treatment (Intercept) 4.033e+00 2.0082238
##  variable:Food_type         (Intercept) 7.592e+00 2.7553900
## Number of obs: 1320, groups:  
## variable, 44; variable:Microbe_treatment, 132; variable:Food_type, 88
## 
## Conditional model:
##                                           Estimate Std. Error z value Pr(>|z|)
## (Intercept)                               -6.88584    0.57879 -11.897  < 2e-16
## Microbe_treatmentMEM                       0.34474    0.48995   0.704  0.48167
## Microbe_treatmentMMO                      -1.70112    0.56011  -3.037  0.00239
## Food_typeLarval food                       0.05614    0.68107   0.082  0.93431
## Microbe_treatmentMEM:Food_typeLarval food -0.47751    0.28648  -1.667  0.09555
## Microbe_treatmentMMO:Food_typeLarval food  0.04079    0.33684   0.121  0.90361
##                                              
## (Intercept)                               ***
## Microbe_treatmentMEM                         
## Microbe_treatmentMMO                      ** 
## Food_typeLarval food                         
## Microbe_treatmentMEM:Food_typeLarval food .  
## Microbe_treatmentMMO:Food_typeLarval food    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Dispersion model:
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)           5.66150    0.58618   9.658  < 2e-16 ***
## variableASV0003       2.32708    0.63327   3.675 0.000238 ***
## variableASV0005       1.09430    0.85645   1.278 0.201352    
## variableASV0006       2.73578    0.66927   4.088 4.36e-05 ***
## variableASV0007       1.53363    0.66346   2.312 0.020803 *  
## variableASV0009       1.95568    0.69320   2.821 0.004784 ** 
## variableASV0010       0.86886    0.73524   1.182 0.237309    
## variableASV0011      -0.13682    0.75431  -0.181 0.856065    
## variableASV0012       1.79164    1.01061   1.773 0.076255 .  
## variableASV0013       2.81004    0.76910   3.654 0.000259 ***
## variableASV0016       1.53770    0.86280   1.782 0.074714 .  
## variableASV0019       2.77594    0.71214   3.898 9.70e-05 ***
## variableASV0024       1.27234    0.70603   1.802 0.071526 .  
## variableASV0027       1.88812    0.73373   2.573 0.010073 *  
## variableASV0032       1.93115    0.79952   2.415 0.015718 *  
## variableASV0033       2.63619    0.72605   3.631 0.000282 ***
## variableASV0043       4.33203    0.88658   4.886 1.03e-06 ***
## variableASV0053       1.06864    0.73318   1.458 0.144966    
## variableASV0054       1.78701    0.92632   1.929 0.053713 .  
## variableASV0055       2.28938    0.78728   2.908 0.003638 ** 
## variableASV0058       2.83787    0.84856   3.344 0.000825 ***
## variableASV0059       0.57687    0.71151   0.811 0.417494    
## variableASV0060       3.00739    0.83668   3.594 0.000325 ***
## variableASV0061       1.60711    0.81783   1.965 0.049405 *  
## variableASV0066      -0.48146    0.74257  -0.648 0.516746    
## variableASV0070       1.35318    0.94334   1.434 0.151441    
## variableASV0079       0.66061    0.68510   0.964 0.334922    
## variableASV0083       0.51392    0.68511   0.750 0.453175    
## variableASV0085       1.19338    1.03223   1.156 0.247630    
## variableASV0087       1.69651    0.76889   2.206 0.027353 *  
## variableASV0099       1.64988    0.82800   1.993 0.046305 *  
## variableASV0100       0.99953    0.75172   1.330 0.183631    
## variableASV0102      -0.41224    0.86194  -0.478 0.632457    
## variableASV0125       2.14935    0.88169   2.438 0.014778 *  
## variableASV0130       1.17259    0.73248   1.601 0.109411    
## variableASV0134       2.21104    1.01097   2.187 0.028740 *  
## variableASV0138       0.30410    0.88967   0.342 0.732493    
## variableASV0149      -1.22305    0.71326  -1.715 0.086397 .  
## variableASV0150       0.04696    0.81330   0.058 0.953952    
## variableASV0164      -0.19801    0.73680  -0.269 0.788125    
## variableASV0218       0.07840    0.94312   0.083 0.933754    
## variableASV0262       0.71369    1.03328   0.691 0.489751    
## variableASV0289      -0.25622    0.86861  -0.295 0.768007    
## variableASV0311       0.13310    0.88717   0.150 0.880744    
## Food_typeLarval food  0.84689    0.27662   3.062 0.002202 ** 
## Microbe_treatmentMEM -0.74771    0.23252  -3.216 0.001302 ** 
## Microbe_treatmentMMO  0.05841    0.31114   0.188 0.851087    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(larv.mod.nb1,larv.mod.nb1.nofm) #sig **
## Data: df.larv.long
## Models:
## larv.mod.nb1.nofm: value ~ Microbe_treatment + Food_type + (1 | variable) + (1 | , zi=~0, disp=~variable + Food_type + Microbe_treatment
## larv.mod.nb1.nofm:     variable:Microbe_treatment) + (1 | variable:Food_type) + , zi=~0, disp=~variable + Food_type + Microbe_treatment
## larv.mod.nb1.nofm:     Microbe_treatment:Food_type + offset(log(sum)), zi=~0, disp=~variable + Food_type + Microbe_treatment
## larv.mod.nb1: value ~ offset(log(sum)) + Microbe_treatment * Food_type + (1 | , zi=~0, disp=~variable + Food_type + Microbe_treatment
## larv.mod.nb1:     variable) + (1 | variable:Microbe_treatment) + (1 | variable:Food_type) + , zi=~0, disp=~variable + Food_type + Microbe_treatment
## larv.mod.nb1:     (1 | variable:Microbe_treatment:Food_type), zi=~0, disp=~variable + Food_type + Microbe_treatment
##                   Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)   
## larv.mod.nb1.nofm 56 9380.4 9670.8 -4634.2   9268.4                            
## larv.mod.nb1      57 9374.9 9670.5 -4630.5   9260.9 7.5229      1   0.006092 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# ##removing just food
# larv.mod.nb1.nofood <- update(larv.mod.nb1,.~.-(1|variable:Food_type))
# anova(larv.mod.nb1,larv.mod.nb1.nofood) #super sig***
# 
# ##removing just microbes
# larv.mod.nb1.nomicr <- update(larv.mod.nb1,.~.-(1|variable:Microbe_treatment))
# anova(larv.mod.nb1,larv.mod.nb1.nomicr) #super sig***

##checking out nbinom2
larv.mod.nb2.nofm <- update(larv.mod.nb2,.~.-(1|variable:Microbe_treatment:Food_type))
anova(larv.mod.nb2,larv.mod.nb2.nofm)
## Data: df.larv.long
## Models:
## larv.mod.nb2.nofm: value ~ Microbe_treatment + Food_type + (1 | variable) + (1 | , zi=~0, disp=~variable + Food_type + Microbe_treatment
## larv.mod.nb2.nofm:     variable:Microbe_treatment) + (1 | variable:Food_type) + , zi=~0, disp=~variable + Food_type + Microbe_treatment
## larv.mod.nb2.nofm:     Microbe_treatment:Food_type + offset(log(sum)), zi=~0, disp=~variable + Food_type + Microbe_treatment
## larv.mod.nb2: value ~ offset(log(sum)) + Microbe_treatment * Food_type + (1 | , zi=~0, disp=~variable + Food_type + Microbe_treatment
## larv.mod.nb2:     variable) + (1 | variable:Microbe_treatment) + (1 | variable:Food_type) + , zi=~0, disp=~variable + Food_type + Microbe_treatment
## larv.mod.nb2:     (1 | variable:Microbe_treatment:Food_type), zi=~0, disp=~variable + Food_type + Microbe_treatment
##                   Df    AIC   BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
## larv.mod.nb2.nofm 56 9763.9 10054 -4825.9   9651.9                             
## larv.mod.nb2      57 9752.6 10048 -4819.3   9638.6 13.281      1  0.0002681 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##sig**

##checking out binomial
larv.mod.bin.nofm <- update(larv.mod.bin,.~.-(1|variable:Microbe_treatment:Food_type))
summary(larv.mod.bin.nofm)
##  Family: binomial  ( logit )
## Formula:          
## cbind(value, sum - value) ~ Microbe_treatment + Food_type + (1 |  
##     variable) + (1 | variable:Microbe_treatment) + (1 | variable:Food_type) +  
##     (1 | rowid) + Microbe_treatment:Food_type
## Data: df.larv.long
## 
##      AIC      BIC   logLik deviance df.resid 
##   9979.2  10031.1  -4979.6   9959.2     1310 
## 
## Random effects:
## 
## Conditional model:
##  Groups                     Name        Variance  Std.Dev. 
##  variable                   (Intercept) 9.084e-07 0.0009531
##  variable:Microbe_treatment (Intercept) 3.206e+01 5.6618257
##  variable:Food_type         (Intercept) 4.579e+01 6.7668655
##  rowid                      (Intercept) 1.071e+01 3.2732847
## Number of obs: 1320, groups:  
## variable, 44; variable:Microbe_treatment, 132; variable:Food_type, 88; rowid, 1320
## 
## Conditional model:
##                                           Estimate Std. Error z value Pr(>|z|)
## (Intercept)                               -13.4410     1.4590  -9.213  < 2e-16
## Microbe_treatmentMEM                        3.0963     1.3187   2.348 0.018878
## Microbe_treatmentMMO                       -5.6255     1.4865  -3.784 0.000154
## Food_typeLarval food                       -0.5707     1.6528  -0.345 0.729870
## Microbe_treatmentMEM:Food_typeLarval food  -2.2886     0.7195  -3.181 0.001470
## Microbe_treatmentMMO:Food_typeLarval food  -0.9151     0.9024  -1.014 0.310532
##                                              
## (Intercept)                               ***
## Microbe_treatmentMEM                      *  
## Microbe_treatmentMMO                      ***
## Food_typeLarval food                         
## Microbe_treatmentMEM:Food_typeLarval food ** 
## Microbe_treatmentMMO:Food_typeLarval food    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(larv.mod.bin,larv.mod.bin.nofm) #sig***
## Data: df.larv.long
## Models:
## larv.mod.bin.nofm: cbind(value, sum - value) ~ Microbe_treatment + Food_type + (1 | , zi=~0, disp=~1
## larv.mod.bin.nofm:     variable) + (1 | variable:Microbe_treatment) + (1 | variable:Food_type) + , zi=~0, disp=~1
## larv.mod.bin.nofm:     (1 | rowid) + Microbe_treatment:Food_type, zi=~0, disp=~1
## larv.mod.bin: cbind(value, sum - value) ~ Microbe_treatment * Food_type + (1 | , zi=~0, disp=~1
## larv.mod.bin:     variable) + (1 | variable:Microbe_treatment) + (1 | variable:Food_type) + , zi=~0, disp=~1
## larv.mod.bin:     (1 | variable:Microbe_treatment:Food_type) + (1 | rowid), zi=~0, disp=~1
##                   Df    AIC   BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
## larv.mod.bin.nofm 10 9979.2 10031 -4979.6   9959.2                             
## larv.mod.bin      11 9953.6 10011 -4965.8   9931.6 27.627      1  1.471e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
###save conditional modes
larv.mod.ranef <- as.data.frame(ranef(larv.mod.nb1))
#write.csv(larv.mod.ranef, "larv.mod.ranef.csv")

# library("ggstats")
# ggcoef_multicomponents(larv.mod.nb1)

3.2.1 Nbinom1 vs. nbinom2 things

##note: some errors/warnings if the treatments aren't factors for some reason
mvtab <- ddply(df.larv.long,
               .(variable:Microbe_treatment:Food_type),
               summarise,
               callmean=mean(value),
               callvar=var(value))

q1 <- qplot(callmean,callvar,data=mvtab)
print(q1+
        ## linear (quasi-Poisson/NB1) fit
        geom_smooth(method="lm",formula=y~x-1,colour="pink")+
        ## smooth (loess)
        geom_smooth(colour="red")+
        ## semi-quadratic (NB2/LNP)
        geom_smooth(method="lm",formula=y~I(x^2)+offset(x)-1,colour="purple")+
        ## Poisson (v=m)
        geom_abline(a=0,b=1,lty=2))#+
## Warning in geom_abline(a = 0, b = 1, lty = 2): Ignoring unknown parameters: `a`
## and `b`
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

        #ylim(0,1.75e8)+
        #xlim(0,20000)

3.3 Water - all treatments

str(df.all.water.long)
## 'data.frame':    3450 obs. of  8 variables:
##  $ sum              : num  101895 104735 144484 122710 89143 ...
##  $ Mesocosm_id      : chr  "B_MMO.4" "B_MMO.5" "B_EMO.1" "B_EMO.2" ...
##  $ Microbe_treatment: chr  "MMO" "MMO" "EMO" "EMO" ...
##  $ Food_type        : chr  "Bamboo" "Bamboo" "Bamboo" "Bamboo" ...
##  $ Biomass_day40    : num  2.612 0.468 3.949 0 0.621 ...
##  $ variable         : Factor w/ 69 levels "ASV0001","ASV0002",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ value            : num  0 0 0 0 0 ...
##  $ rowid            : int  1 2 3 4 5 6 7 8 9 10 ...
df.all.water.long$rowid <- as.factor(df.all.water.long$rowid)
df.all.water.long$Microbe_treatment <- as.factor(df.all.water.long$Microbe_treatment)
df.all.water.long$Food_type <- as.factor(df.all.water.long$Food_type)
#df.all.water.long$food.microbes <- as.factor(df.all.water.long$food.microbes)
df.all.water.long$Mesocosm_id <- as.factor(df.all.water.long$Mesocosm_id)
str(df.all.water.long)
## 'data.frame':    3450 obs. of  8 variables:
##  $ sum              : num  101895 104735 144484 122710 89143 ...
##  $ Mesocosm_id      : Factor w/ 50 levels "B_ALL.1","B_ALL.2",..: 24 25 11 12 13 14 15 6 7 8 ...
##  $ Microbe_treatment: Factor w/ 5 levels "ALL","ECO","EMO",..: 5 5 3 3 3 3 3 2 2 2 ...
##  $ Food_type        : Factor w/ 2 levels "Bamboo","Larval food": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Biomass_day40    : num  2.612 0.468 3.949 0 0.621 ...
##  $ variable         : Factor w/ 69 levels "ASV0001","ASV0002",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ value            : num  0 0 0 0 0 ...
##  $ rowid            : Factor w/ 3450 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
##full model
water.mod.nb1 <- glmmTMB(value~offset(log(sum))+Microbe_treatment*Food_type+(1|variable)+(1|variable:Microbe_treatment)+(1|variable:Food_type)+(1|variable:Microbe_treatment:Food_type),family=nbinom1,dispformula=~variable+Microbe_treatment+Food_type,data=df.all.water.long)
summary(water.mod.nb1)
##  Family: nbinom1  ( log )
## Formula:          
## value ~ offset(log(sum)) + Microbe_treatment * Food_type + (1 |  
##     variable) + (1 | variable:Microbe_treatment) + (1 | variable:Food_type) +  
##     (1 | variable:Microbe_treatment:Food_type)
## Dispersion:             ~variable + Microbe_treatment + Food_type
## Data: df.all.water.long
## 
##      AIC      BIC   logLik deviance df.resid 
##  22998.2  23539.0 -11411.1  22822.2     3362 
## 
## Random effects:
## 
## Conditional model:
##  Groups                               Name        Variance Std.Dev.
##  variable                             (Intercept) 0.2286   0.4781  
##  variable:Microbe_treatment           (Intercept) 6.3695   2.5238  
##  variable:Food_type                   (Intercept) 7.8079   2.7943  
##  variable:Microbe_treatment:Food_type (Intercept) 0.2461   0.4961  
## Number of obs: 3450, groups:  
## variable, 69; variable:Microbe_treatment, 345; variable:Food_type, 138; variable:Microbe_treatment:Food_type, 690
## 
## Conditional model:
##                                           Estimate Std. Error z value Pr(>|z|)
## (Intercept)                               -6.44631    0.47687 -13.518  < 2e-16
## Microbe_treatmentECO                      -7.58589    0.80260  -9.452  < 2e-16
## Microbe_treatmentEMO                       0.08668    0.45898   0.189 0.850211
## Microbe_treatmentMEM                      -0.12823    0.45788  -0.280 0.779435
## Microbe_treatmentMMO                      -2.51900    0.57718  -4.364 1.28e-05
## Food_typeLarval food                      -2.66398    0.54486  -4.889 1.01e-06
## Microbe_treatmentECO:Food_typeLarval food -1.34724    0.58473  -2.304 0.021220
## Microbe_treatmentEMO:Food_typeLarval food  0.56380    0.23769   2.372 0.017691
## Microbe_treatmentMEM:Food_typeLarval food  0.90291    0.23485   3.845 0.000121
## Microbe_treatmentMMO:Food_typeLarval food  0.54351    0.37460   1.451 0.146814
##                                              
## (Intercept)                               ***
## Microbe_treatmentECO                      ***
## Microbe_treatmentEMO                         
## Microbe_treatmentMEM                         
## Microbe_treatmentMMO                      ***
## Food_typeLarval food                      ***
## Microbe_treatmentECO:Food_typeLarval food *  
## Microbe_treatmentEMO:Food_typeLarval food *  
## Microbe_treatmentMEM:Food_typeLarval food ***
## Microbe_treatmentMMO:Food_typeLarval food    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Dispersion model:
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)           6.37016    0.42392  15.027  < 2e-16 ***
## variableASV0002       0.77344    0.46332   1.669 0.095052 .  
## variableASV0003      -0.92224    0.55333  -1.667 0.095571 .  
## variableASV0004       2.63337    0.51177   5.146 2.67e-07 ***
## variableASV0005       2.25959    0.50658   4.461 8.18e-06 ***
## variableASV0007       0.97636    0.50479   1.934 0.053090 .  
## variableASV0008       1.31909    0.51861   2.544 0.010974 *  
## variableASV0009       0.21037    0.48832   0.431 0.666613    
## variableASV0010       1.43391    0.50300   2.851 0.004362 ** 
## variableASV0011       1.88834    0.51543   3.664 0.000249 ***
## variableASV0012       1.49855    0.50859   2.947 0.003214 ** 
## variableASV0015       2.47173    0.55981   4.415 1.01e-05 ***
## variableASV0016       0.08453    0.51955   0.163 0.870757    
## variableASV0017       0.09542    0.58335   0.164 0.870062    
## variableASV0018      -0.11272    0.89425  -0.126 0.899692    
## variableASV0020       1.87992    0.53161   3.536 0.000406 ***
## variableASV0021       1.52117    0.76308   1.993 0.046210 *  
## variableASV0022       1.07533    0.51278   2.097 0.035988 *  
## variableASV0023       1.46196    0.61914   2.361 0.018212 *  
## variableASV0024       1.40504    0.52409   2.681 0.007342 ** 
## variableASV0025      -0.21316    0.60895  -0.350 0.726302    
## variableASV0027      -0.73196    0.59584  -1.228 0.219277    
## variableASV0033       1.39685    0.96018   1.455 0.145731    
## variableASV0036       3.75485    0.67649   5.550 2.85e-08 ***
## variableASV0037       3.54829    0.86838   4.086 4.39e-05 ***
## variableASV0039       0.73133    0.58984   1.240 0.215019    
## variableASV0040       0.47978    0.69209   0.693 0.488166    
## variableASV0042       2.09878    0.54754   3.833 0.000127 ***
## variableASV0044       1.26503    0.59555   2.124 0.033659 *  
## variableASV0051      -0.50545    0.50943  -0.992 0.321103    
## variableASV0054       0.49609    0.71855   0.690 0.489947    
## variableASV0057       0.14685    0.61464   0.239 0.811161    
## variableASV0059       0.68076    0.58675   1.160 0.245959    
## variableASV0063       0.26277    0.54803   0.479 0.631592    
## variableASV0065      -0.30281    0.61917  -0.489 0.624801    
## variableASV0066       1.51543    0.72280   2.097 0.036030 *  
## variableASV0069      -1.69629    0.62018  -2.735 0.006235 ** 
## variableASV0072      -0.43166    0.59680  -0.723 0.469506    
## variableASV0073       1.70443    0.67058   2.542 0.011031 *  
## variableASV0079      -0.49152    0.78820  -0.624 0.532894    
## variableASV0080       1.69708    0.69186   2.453 0.014170 *  
## variableASV0082       1.15666    0.63788   1.813 0.069785 .  
## variableASV0083      -1.25936    0.81448  -1.546 0.122054    
## variableASV0084       0.48664    0.62762   0.775 0.438117    
## variableASV0085      -0.18806    0.53040  -0.355 0.722916    
## variableASV0089       1.01164    0.65379   1.547 0.121782    
## variableASV0091       1.06317    0.80740   1.317 0.187911    
## variableASV0093      -1.40509    0.60456  -2.324 0.020117 *  
## variableASV0094      -0.42514    0.60416  -0.704 0.481628    
## variableASV0095       0.50560    0.63661   0.794 0.427074    
## variableASV0101       1.64515    0.80063   2.055 0.039897 *  
## variableASV0102       0.34004    0.68938   0.493 0.621834    
## variableASV0105       1.59824    1.07869   1.482 0.138431    
## variableASV0111      -0.95686    0.57235  -1.672 0.094560 .  
## variableASV0112       1.24346    0.76621   1.623 0.104617    
## variableASV0115       0.53222    0.64288   0.828 0.407742    
## variableASV0116       0.03398    0.57100   0.060 0.952550    
## variableASV0119      -0.03823    0.62799  -0.061 0.951459    
## variableASV0123       0.21105    0.65198   0.324 0.746158    
## variableASV0128       1.26038    0.83501   1.509 0.131194    
## variableASV0139       1.53228    0.82633   1.854 0.063694 .  
## variableASV0145      -1.86479    0.94412  -1.975 0.048250 *  
## variableASV0146       0.66931    0.80331   0.833 0.404737    
## variableASV0152       1.39689    1.10004   1.270 0.204138    
## variableASV0155      -1.17553    0.71711  -1.639 0.101157    
## variableASV0167       0.31581    0.69560   0.454 0.649816    
## variableASV0175       1.43254    1.10308   1.299 0.194058    
## variableASV0180      -1.07281    0.62117  -1.727 0.084152 .  
## variableASV0184       1.09999    0.98249   1.120 0.262889    
## Microbe_treatmentECO -1.33148    0.52429  -2.540 0.011098 *  
## Microbe_treatmentEMO -0.59136    0.14822  -3.990 6.61e-05 ***
## Microbe_treatmentMEM -0.66870    0.14362  -4.656 3.22e-06 ***
## Microbe_treatmentMMO  1.35596    0.24013   5.647 1.63e-08 ***
## Food_typeLarval food -0.29733    0.13983  -2.126 0.033470 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
water.mod.nb2 <- glmmTMB(value~offset(log(sum))+Microbe_treatment*Food_type+(1|variable)+(1|variable:Microbe_treatment)+(1|variable:Food_type)+(1|variable:Microbe_treatment:Food_type),family=nbinom2,dispformula=~variable+Microbe_treatment+Food_type,data=df.all.water.long)
summary(water.mod.nb2)
##  Family: nbinom2  ( log )
## Formula:          
## value ~ offset(log(sum)) + Microbe_treatment * Food_type + (1 |  
##     variable) + (1 | variable:Microbe_treatment) + (1 | variable:Food_type) +  
##     (1 | variable:Microbe_treatment:Food_type)
## Dispersion:             ~variable + Microbe_treatment + Food_type
## Data: df.all.water.long
## 
##      AIC      BIC   logLik deviance df.resid 
##  23799.3  24340.2 -11811.7  23623.3     3362 
## 
## Random effects:
## 
## Conditional model:
##  Groups                               Name        Variance  Std.Dev.
##  variable                             (Intercept) 3.444e-06 0.001856
##  variable:Microbe_treatment           (Intercept) 1.963e+01 4.431121
##  variable:Food_type                   (Intercept) 2.196e+01 4.686453
##  variable:Microbe_treatment:Food_type (Intercept) 3.387e-01 0.582020
## Number of obs: 3450, groups:  
## variable, 69; variable:Microbe_treatment, 345; variable:Food_type, 138; variable:Microbe_treatment:Food_type, 690
## 
## Conditional model:
##                                            Estimate Std. Error z value Pr(>|z|)
## (Intercept)                                -6.70990    0.80867  -8.297  < 2e-16
## Microbe_treatmentECO                      -13.46838    0.98961 -13.610  < 2e-16
## Microbe_treatmentEMO                       -0.07208    0.80083  -0.090 0.928285
## Microbe_treatmentMEM                       -0.25176    0.79989  -0.315 0.752957
## Microbe_treatmentMMO                       -6.68505    1.22188  -5.471 4.47e-08
## Food_typeLarval food                       -4.75081    0.91058  -5.217 1.82e-07
## Microbe_treatmentECO:Food_typeLarval food  -1.41242    0.78679  -1.795 0.072627
## Microbe_treatmentEMO:Food_typeLarval food   1.24419    0.37240   3.341 0.000835
## Microbe_treatmentMEM:Food_typeLarval food   1.51652    0.37018   4.097 4.19e-05
## Microbe_treatmentMMO:Food_typeLarval food   1.19207    0.60428   1.973 0.048528
##                                              
## (Intercept)                               ***
## Microbe_treatmentECO                      ***
## Microbe_treatmentEMO                         
## Microbe_treatmentMEM                         
## Microbe_treatmentMMO                      ***
## Food_typeLarval food                      ***
## Microbe_treatmentECO:Food_typeLarval food .  
## Microbe_treatmentEMO:Food_typeLarval food ***
## Microbe_treatmentMEM:Food_typeLarval food ***
## Microbe_treatmentMMO:Food_typeLarval food *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Dispersion model:
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)          -1.59171    0.35208  -4.521 6.16e-06 ***
## variableASV0002       2.98368    0.40931   7.289 3.11e-13 ***
## variableASV0003       2.29908    0.45337   5.071 3.95e-07 ***
## variableASV0004       2.69638    0.42981   6.273 3.53e-10 ***
## variableASV0005       3.09013    0.43441   7.113 1.13e-12 ***
## variableASV0007       1.94350    0.41292   4.707 2.52e-06 ***
## variableASV0008       1.44625    0.42608   3.394 0.000688 ***
## variableASV0009       3.61752    0.40976   8.828  < 2e-16 ***
## variableASV0010       2.95017    0.43147   6.837 8.06e-12 ***
## variableASV0011       1.53679    0.42889   3.583 0.000339 ***
## variableASV0012       1.67829    0.41171   4.076 4.57e-05 ***
## variableASV0015       0.32703    0.41830   0.782 0.434336    
## variableASV0016       4.41497    0.43980  10.039  < 2e-16 ***
## variableASV0017       1.59587    0.49810   3.204 0.001356 ** 
## variableASV0018       0.20353    0.61338   0.332 0.740028    
## variableASV0020       1.26221    0.52218   2.417 0.015640 *  
## variableASV0021      -0.19304    0.53055  -0.364 0.715974    
## variableASV0022       1.59617    0.42210   3.781 0.000156 ***
## variableASV0023       0.94622    0.51645   1.832 0.066928 .  
## variableASV0024       1.55626    0.42474   3.664 0.000248 ***
## variableASV0025       3.07164    0.52663   5.833 5.46e-09 ***
## variableASV0027       2.27425    0.54024   4.210 2.56e-05 ***
## variableASV0033      -1.70585    0.55628  -3.067 0.002165 ** 
## variableASV0036      -0.89001    0.44272  -2.010 0.044398 *  
## variableASV0037      -2.02357    0.49563  -4.083 4.45e-05 ***
## variableASV0039       0.16350    0.54971   0.297 0.766140    
## variableASV0040       0.24090    0.52449   0.459 0.646016    
## variableASV0042       0.70825    0.41900   1.690 0.090962 .  
## variableASV0044       1.97892    0.50990   3.881 0.000104 ***
## variableASV0051       1.34829    0.42171   3.197 0.001388 ** 
## variableASV0054      -0.79006    0.49229  -1.605 0.108520    
## variableASV0057       0.08184    0.48140   0.170 0.865012    
## variableASV0059      -0.03230    0.45870  -0.070 0.943862    
## variableASV0063       0.51779    0.44137   1.173 0.240743    
## variableASV0065       0.62284    0.49265   1.264 0.206137    
## variableASV0066      -0.77537    0.49398  -1.570 0.116499    
## variableASV0069       3.97236    0.53702   7.397 1.39e-13 ***
## variableASV0072       2.89809    0.52554   5.515 3.50e-08 ***
## variableASV0073      -0.81592    0.47620  -1.713 0.086637 .  
## variableASV0079       0.29782    0.55622   0.535 0.592347    
## variableASV0080      -0.84527    0.47211  -1.790 0.073391 .  
## variableASV0082      -0.52171    0.45860  -1.138 0.255277    
## variableASV0083       1.51843    0.59704   2.543 0.010982 *  
## variableASV0084       0.40152    0.47282   0.849 0.395763    
## variableASV0085       0.94584    0.43706   2.164 0.030455 *  
## variableASV0089       0.58796    0.51865   1.134 0.256947    
## variableASV0091      -0.82137    0.56218  -1.461 0.144006    
## variableASV0093       2.76554    0.52528   5.265 1.40e-07 ***
## variableASV0094       3.42511    0.51985   6.589 4.44e-11 ***
## variableASV0095       0.90323    0.51196   1.764 0.077688 .  
## variableASV0101      -1.21368    0.48963  -2.479 0.013184 *  
## variableASV0102       0.73940    0.51083   1.447 0.147773    
## variableASV0105      -1.04665    0.74319  -1.408 0.159036    
## variableASV0111       0.08368    0.46172   0.181 0.856185    
## variableASV0112      -0.23926    0.52980  -0.452 0.651548    
## variableASV0115      -0.48014    0.49760  -0.965 0.334588    
## variableASV0116      -0.11382    0.43731  -0.260 0.794664    
## variableASV0119       1.15229    0.52336   2.202 0.027686 *  
## variableASV0123       0.64408    0.52161   1.235 0.216909    
## variableASV0128      -1.33402    0.52480  -2.542 0.011024 *  
## variableASV0139      -0.57802    0.54388  -1.063 0.287886    
## variableASV0145       1.77656    0.63592   2.794 0.005211 ** 
## variableASV0146      -0.72418    0.54123  -1.338 0.180886    
## variableASV0152      -1.40699    0.60350  -2.331 0.019733 *  
## variableASV0155       3.57915    0.59920   5.973 2.33e-09 ***
## variableASV0167       0.36274    0.52882   0.686 0.492749    
## variableASV0175      -1.40188    0.60345  -2.323 0.020173 *  
## variableASV0180       1.72267    0.54713   3.149 0.001641 ** 
## variableASV0184      -1.83812    0.55782  -3.295 0.000984 ***
## Microbe_treatmentECO  3.59384    0.42409   8.474  < 2e-16 ***
## Microbe_treatmentEMO  0.60138    0.12212   4.924 8.46e-07 ***
## Microbe_treatmentMEM  0.49914    0.11656   4.282 1.85e-05 ***
## Microbe_treatmentMMO -1.04862    0.23045  -4.550 5.36e-06 ***
## Food_typeLarval food -0.50118    0.12017  -4.171 3.04e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
AICtab(water.mod.nb1,water.mod.nb2)
##               dAIC  df
## water.mod.nb1   0.0 88
## water.mod.nb2 801.1 88
water.mod.nb1.resid <- simulateResiduals(fittedModel = water.mod.nb1, plot = T)

water.mod.nb2.resid <- simulateResiduals(fittedModel = water.mod.nb2, plot = T)
## DHARMa:testOutliers with type = binomial may have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details

##dispersal formula checks
water.mod.nb1.disnovar <- glmmTMB(value~offset(log(sum))+Microbe_treatment*Food_type+(1|variable)+(1|variable:Microbe_treatment)+(1|variable:Food_type)+(1|variable:Microbe_treatment:Food_type),family=nbinom1,dispformula=~Microbe_treatment+Food_type,data=df.all.water.long)

anova(water.mod.nb1.disnovar,water.mod.nb1) #sig***
## Data: df.all.water.long
## Models:
## water.mod.nb1.disnovar: value ~ offset(log(sum)) + Microbe_treatment * Food_type + (1 | , zi=~0, disp=~Microbe_treatment + Food_type
## water.mod.nb1.disnovar:     variable) + (1 | variable:Microbe_treatment) + (1 | variable:Food_type) + , zi=~0, disp=~variable + Microbe_treatment + Food_type
## water.mod.nb1.disnovar:     (1 | variable:Microbe_treatment:Food_type), zi=~0, disp=~Microbe_treatment + Food_type
## water.mod.nb1: value ~ offset(log(sum)) + Microbe_treatment * Food_type + (1 | , zi=~0, disp=~variable + Microbe_treatment + Food_type
## water.mod.nb1:     variable) + (1 | variable:Microbe_treatment) + (1 | variable:Food_type) + , zi=~0, disp=~Microbe_treatment + Food_type
## water.mod.nb1:     (1 | variable:Microbe_treatment:Food_type), zi=~0, disp=~variable + Microbe_treatment + Food_type
##                        Df   AIC   BIC logLik deviance  Chisq Chi Df Pr(>Chisq)
## water.mod.nb1.disnovar 20 23320 23443 -11640    23280                         
## water.mod.nb1          88 22998 23539 -11411    22822 458.29     68  < 2.2e-16
##                           
## water.mod.nb1.disnovar    
## water.mod.nb1          ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##removing 3 way interaction
water.mod.nb1.nofm <- update(water.mod.nb1,.~.-(1|variable:Microbe_treatment:Food_type))
anova(water.mod.nb1,water.mod.nb1.nofm)
## Data: df.all.water.long
## Models:
## water.mod.nb1.nofm: value ~ Microbe_treatment + Food_type + (1 | variable) + (1 | , zi=~0, disp=~variable + Microbe_treatment + Food_type
## water.mod.nb1.nofm:     variable:Microbe_treatment) + (1 | variable:Food_type) + , zi=~0, disp=~variable + Microbe_treatment + Food_type
## water.mod.nb1.nofm:     Microbe_treatment:Food_type + offset(log(sum)), zi=~0, disp=~variable + Microbe_treatment + Food_type
## water.mod.nb1: value ~ offset(log(sum)) + Microbe_treatment * Food_type + (1 | , zi=~0, disp=~variable + Microbe_treatment + Food_type
## water.mod.nb1:     variable) + (1 | variable:Microbe_treatment) + (1 | variable:Food_type) + , zi=~0, disp=~variable + Microbe_treatment + Food_type
## water.mod.nb1:     (1 | variable:Microbe_treatment:Food_type), zi=~0, disp=~variable + Microbe_treatment + Food_type
##                    Df   AIC   BIC logLik deviance  Chisq Chi Df Pr(>Chisq)    
## water.mod.nb1.nofm 87 23065 23600 -11446    22891                             
## water.mod.nb1      88 22998 23539 -11411    22822 69.225      1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##super sig 0.000242

water.mod.nb2.nofm <- update(water.mod.nb2,.~.-(1|variable:Microbe_treatment:Food_type))
anova(water.mod.nb2,water.mod.nb2.nofm)
## Data: df.all.water.long
## Models:
## water.mod.nb2.nofm: value ~ Microbe_treatment + Food_type + (1 | variable) + (1 | , zi=~0, disp=~variable + Microbe_treatment + Food_type
## water.mod.nb2.nofm:     variable:Microbe_treatment) + (1 | variable:Food_type) + , zi=~0, disp=~variable + Microbe_treatment + Food_type
## water.mod.nb2.nofm:     Microbe_treatment:Food_type + offset(log(sum)), zi=~0, disp=~variable + Microbe_treatment + Food_type
## water.mod.nb2: value ~ offset(log(sum)) + Microbe_treatment * Food_type + (1 | , zi=~0, disp=~variable + Microbe_treatment + Food_type
## water.mod.nb2:     variable) + (1 | variable:Microbe_treatment) + (1 | variable:Food_type) + , zi=~0, disp=~variable + Microbe_treatment + Food_type
## water.mod.nb2:     (1 | variable:Microbe_treatment:Food_type), zi=~0, disp=~variable + Microbe_treatment + Food_type
##                    Df   AIC   BIC logLik deviance  Chisq Chi Df Pr(>Chisq)    
## water.mod.nb2.nofm 87 23825 24360 -11826    23651                             
## water.mod.nb2      88 23799 24340 -11812    23623 27.894      1  1.281e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##also super sig

###save conditional modes
water.mod.ranef <- as.data.frame(ranef(water.mod.nb1))
#write.csv(water.mod.ranef, "water.mod.ranef.all.csv")
#saveRDS(water.mod.nb1,file="water.mod.nb1.all.rds")

3.3.1 Nbinom1 vs. nbinom2 things

Following this wonderful resource, on page 16

Ben Bolker, Mollie Brooks, Beth Gardner, Cleridy Lennert, Mihoko Minami, October 23, 2012, Owls example: a zero-inflated, generalized linear mixed model for count data. https://groups.nceas.ucsb.edu/non-linear-modeling/projects/owls/WRITEUP/owls.pdf

mvtab <- ddply(df.all.water.long,
               .(variable:Microbe_treatment:Food_type),
               summarise,
               callmean=mean(value),
               callvar=var(value))

q1 <- qplot(callmean,callvar,data=mvtab)
print(q1+
        ## linear (quasi-Poisson/NB1) fit
        geom_smooth(method="lm",formula=y~x-1,colour="pink")+
        ## smooth (loess)
        geom_smooth(colour="red")+
        ## semi-quadratic (NB2/LNP)
        geom_smooth(method="lm",formula=y~I(x^2)+offset(x)-1,colour="purple")+
        ## Poisson (v=m)
        geom_abline(a=0,b=1,lty=2))
## Warning in geom_abline(a = 0, b = 1, lty = 2): Ignoring unknown parameters: `a`
## and `b`
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

3.4 Larvae - all treatments

str(df.all.larv.long)
## 'data.frame':    2700 obs. of  8 variables:
##  $ sum              : num  63934 18777 13503 111 2064 ...
##  $ Mesocosm_id      : chr  "B_EMO.3" "B_EMO.4" "B_EMO.5" "B_ECO.1" ...
##  $ Microbe_treatment: chr  "EMO" "EMO" "EMO" "ECO" ...
##  $ Food_type        : chr  "Bamboo" "Bamboo" "Bamboo" "Bamboo" ...
##  $ Biomass_day40    : num  0.621 0.308 4.085 0 0 ...
##  $ variable         : Factor w/ 54 levels "ASV0001","ASV0002",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ value            : num  0 0 15 111 1867 ...
##  $ rowid            : int  1 2 3 4 5 6 7 8 9 10 ...
df.all.larv.long$rowid <- as.factor(df.all.larv.long$rowid)
df.all.larv.long$Microbe_treatment <- as.factor(df.all.larv.long$Microbe_treatment)
df.all.larv.long$Food_type <- as.factor(df.all.larv.long$Food_type)
#df.all.larv.long$food.microbes <- as.factor(df.all.larv.long$food.microbes)
df.all.larv.long$Mesocosm_id <- as.factor(df.all.larv.long$Mesocosm_id)
str(df.all.larv.long)
## 'data.frame':    2700 obs. of  8 variables:
##  $ sum              : num  63934 18777 13503 111 2064 ...
##  $ Mesocosm_id      : Factor w/ 50 levels "B_ALL.1","B_ALL.2",..: 13 14 15 6 7 8 9 10 16 17 ...
##  $ Microbe_treatment: Factor w/ 5 levels "ALL","ECO","EMO",..: 3 3 3 2 2 2 2 2 4 4 ...
##  $ Food_type        : Factor w/ 2 levels "Bamboo","Larval food": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Biomass_day40    : num  0.621 0.308 4.085 0 0 ...
##  $ variable         : Factor w/ 54 levels "ASV0001","ASV0002",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ value            : num  0 0 15 111 1867 ...
##  $ rowid            : Factor w/ 2700 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
##binomial
larv.mod.bin <- glmmTMB(cbind(value, sum-value)~Microbe_treatment*Food_type+(1|variable)+(1|variable:Microbe_treatment)+(1|variable:Food_type)+(1|variable:Microbe_treatment:Food_type)+(1|rowid),family=binomial,data=df.all.larv.long)
summary(larv.mod.bin)
##  Family: binomial  ( logit )
## Formula:          
## cbind(value, sum - value) ~ Microbe_treatment * Food_type + (1 |  
##     variable) + (1 | variable:Microbe_treatment) + (1 | variable:Food_type) +  
##     (1 | variable:Microbe_treatment:Food_type) + (1 | rowid)
## Data: df.all.larv.long
## 
##      AIC      BIC   logLik deviance df.resid 
##  15583.1  15671.6  -7776.6  15553.1     2685 
## 
## Random effects:
## 
## Conditional model:
##  Groups                               Name        Variance  Std.Dev. 
##  variable                             (Intercept) 8.950e-07 0.0009461
##  variable:Microbe_treatment           (Intercept) 2.604e+01 5.1032341
##  variable:Food_type                   (Intercept) 3.717e+01 6.0969133
##  variable:Microbe_treatment:Food_type (Intercept) 6.807e+00 2.6090817
##  rowid                                (Intercept) 1.322e+01 3.6365737
## Number of obs: 2700, groups:  
## variable, 54; variable:Microbe_treatment, 270; variable:Food_type, 108; variable:Microbe_treatment:Food_type, 540; rowid, 2700
## 
## Conditional model:
##                                           Estimate Std. Error z value Pr(>|z|)
## (Intercept)                               -11.4884     1.2226  -9.397  < 2e-16
## Microbe_treatmentECO                       -8.9465     1.3737  -6.513 7.37e-11
## Microbe_treatmentEMO                       -2.8262     1.2533  -2.255 0.024137
## Microbe_treatmentMEM                        0.1802     1.2341   0.146 0.883936
## Microbe_treatmentMMO                       -8.2294     1.2300  -6.691 2.22e-11
## Food_typeLarval food                       -5.8031     1.4682  -3.953 7.73e-05
## Microbe_treatmentECO:Food_typeLarval food  -2.3232     1.8710  -1.242 0.214342
## Microbe_treatmentEMO:Food_typeLarval food   3.8919     1.1002   3.538 0.000404
## Microbe_treatmentMEM:Food_typeLarval food   1.6363     1.0968   1.492 0.135748
## Microbe_treatmentMMO:Food_typeLarval food   2.8462     1.2093   2.354 0.018593
##                                              
## (Intercept)                               ***
## Microbe_treatmentECO                      ***
## Microbe_treatmentEMO                      *  
## Microbe_treatmentMEM                         
## Microbe_treatmentMMO                      ***
## Food_typeLarval food                      ***
## Microbe_treatmentECO:Food_typeLarval food    
## Microbe_treatmentEMO:Food_typeLarval food ***
## Microbe_treatmentMEM:Food_typeLarval food    
## Microbe_treatmentMMO:Food_typeLarval food *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##full model
larv.mod.nb1 <- glmmTMB(value~offset(log(sum))+Microbe_treatment*Food_type+(1|variable)+(1|variable:Microbe_treatment)+(1|variable:Food_type)+(1|variable:Microbe_treatment:Food_type),disp=~variable+Food_type+Microbe_treatment,family=nbinom1,data=df.all.larv.long)
summary(larv.mod.nb1)
##  Family: nbinom1  ( log )
## Formula:          
## value ~ offset(log(sum)) + Microbe_treatment * Food_type + (1 |  
##     variable) + (1 | variable:Microbe_treatment) + (1 | variable:Food_type) +  
##     (1 | variable:Microbe_treatment:Food_type)
## Dispersion:             ~variable + Food_type + Microbe_treatment
## Data: df.all.larv.long
## 
##      AIC      BIC   logLik deviance df.resid 
##  14676.1  15106.9  -7265.0  14530.1     2627 
## 
## Random effects:
## 
## Conditional model:
##  Groups                               Name        Variance  Std.Dev. 
##  variable                             (Intercept) 1.860e-07 0.0004313
##  variable:Microbe_treatment           (Intercept) 3.833e+00 1.9577005
##  variable:Food_type                   (Intercept) 6.559e+00 2.5610822
##  variable:Microbe_treatment:Food_type (Intercept) 2.367e-01 0.4864757
## Number of obs: 2700, groups:  
## variable, 54; variable:Microbe_treatment, 270; variable:Food_type, 108; variable:Microbe_treatment:Food_type, 540
## 
## Conditional model:
##                                           Estimate Std. Error z value Pr(>|z|)
## (Intercept)                                -6.4933     0.4805 -13.513  < 2e-16
## Microbe_treatmentECO                       -2.8758     0.9653  -2.979 0.002889
## Microbe_treatmentEMO                       -0.5272     0.4613  -1.143 0.253106
## Microbe_treatmentMEM                       -0.1958     0.4377  -0.447 0.654578
## Microbe_treatmentMMO                       -2.1999     0.5132  -4.287 1.81e-05
## Food_typeLarval food                       -1.1961     0.5968  -2.004 0.045035
## Microbe_treatmentECO:Food_typeLarval food  -2.3653     0.6487  -3.646 0.000266
## Microbe_treatmentEMO:Food_typeLarval food   0.9873     0.3431   2.877 0.004010
## Microbe_treatmentMEM:Food_typeLarval food   0.7123     0.3163   2.252 0.024326
## Microbe_treatmentMMO:Food_typeLarval food   1.1627     0.3863   3.010 0.002611
##                                              
## (Intercept)                               ***
## Microbe_treatmentECO                      ** 
## Microbe_treatmentEMO                         
## Microbe_treatmentMEM                         
## Microbe_treatmentMMO                      ***
## Food_typeLarval food                      *  
## Microbe_treatmentECO:Food_typeLarval food ***
## Microbe_treatmentEMO:Food_typeLarval food ** 
## Microbe_treatmentMEM:Food_typeLarval food *  
## Microbe_treatmentMMO:Food_typeLarval food ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Dispersion model:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)           6.608016   0.480928  13.740  < 2e-16 ***
## variableASV0002       0.125802   0.830243   0.152 0.879562    
## variableASV0003       1.424233   0.541764   2.629 0.008567 ** 
## variableASV0004       1.377207   0.989981   1.391 0.164181    
## variableASV0005       0.568316   0.729597   0.779 0.436012    
## variableASV0006       1.481056   0.556652   2.661 0.007799 ** 
## variableASV0007       0.205723   0.574432   0.358 0.720244    
## variableASV0009       0.805902   0.583140   1.382 0.166971    
## variableASV0010       0.120688   0.643408   0.188 0.851209    
## variableASV0011      -0.744663   0.627092  -1.187 0.235036    
## variableASV0012       0.397982   0.814716   0.488 0.625202    
## variableASV0013       1.423164   0.668655   2.128 0.033304 *  
## variableASV0016       0.710856   0.765013   0.929 0.352781    
## variableASV0019       1.814936   0.676409   2.683 0.007292 ** 
## variableASV0024      -0.120270   0.592496  -0.203 0.839144    
## variableASV0027       0.616208   0.640521   0.962 0.336029    
## variableASV0030      -0.101724   1.076445  -0.095 0.924712    
## variableASV0032       0.452840   0.710405   0.637 0.523839    
## variableASV0033       1.277173   0.603836   2.115 0.034422 *  
## variableASV0043       3.361423   0.847375   3.967 7.28e-05 ***
## variableASV0053       0.366468   0.629382   0.582 0.560387    
## variableASV0054       0.409139   0.820909   0.498 0.618204    
## variableASV0055       1.455563   0.696111   2.091 0.036529 *  
## variableASV0058       1.669642   0.711720   2.346 0.018980 *  
## variableASV0059      -0.524573   0.607181  -0.864 0.387616    
## variableASV0060       2.074881   0.759824   2.731 0.006319 ** 
## variableASV0061       1.279328   0.660277   1.938 0.052677 .  
## variableASV0066      -1.241301   0.686333  -1.809 0.070513 .  
## variableASV0070      -0.060359   0.812150  -0.074 0.940756    
## variableASV0079      -0.709432   0.587382  -1.208 0.227129    
## variableASV0083      -0.604105   0.588541  -1.026 0.304681    
## variableASV0085       0.009405   0.958822   0.010 0.992174    
## variableASV0087       0.545854   0.655877   0.832 0.405268    
## variableASV0099       0.455684   0.725616   0.628 0.530006    
## variableASV0100      -0.289356   0.659121  -0.439 0.660660    
## variableASV0102      -1.381170   0.862920  -1.601 0.109471    
## variableASV0109       0.596347   1.010469   0.590 0.555078    
## variableASV0125       1.346848   0.776236   1.735 0.082723 .  
## variableASV0130       0.076914   0.678606   0.113 0.909760    
## variableASV0134       0.931249   0.904047   1.030 0.302968    
## variableASV0138      -0.864209   0.741853  -1.165 0.244046    
## variableASV0149      -2.035605   0.624026  -3.262 0.001106 ** 
## variableASV0150      -1.340375   0.689884  -1.943 0.052028 .  
## variableASV0152       0.891285   0.907213   0.982 0.325882    
## variableASV0164      -1.121652   0.661599  -1.695 0.090006 .  
## variableASV0210       0.365497   1.072546   0.341 0.733273    
## variableASV0212       0.620222   0.959904   0.646 0.518196    
## variableASV0218      -1.112410   0.845764  -1.315 0.188418    
## variableASV0220       1.082720   1.052458   1.029 0.303595    
## variableASV0225       0.297385   1.068749   0.278 0.780816    
## variableASV0262      -0.372923   1.053183  -0.354 0.723271    
## variableASV0276       0.684188   1.024039   0.668 0.504053    
## variableASV0289      -1.178319   0.801413  -1.470 0.141480    
## variableASV0311      -0.957374   0.776274  -1.233 0.217466    
## Food_typeLarval food  1.045981   0.222090   4.710 2.48e-06 ***
## Microbe_treatmentECO -1.008090   0.807555  -1.248 0.211913    
## Microbe_treatmentEMO -0.100263   0.223110  -0.449 0.653152    
## Microbe_treatmentMEM -0.661644   0.194947  -3.394 0.000689 ***
## Microbe_treatmentMMO  0.208428   0.263089   0.792 0.428224    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
larv.mod.nb2 <- glmmTMB(value~offset(log(sum))+Microbe_treatment*Food_type+(1|variable)+(1|variable:Microbe_treatment)+(1|variable:Food_type)+(1|variable:Microbe_treatment:Food_type),disp=~variable+Food_type+Microbe_treatment,family=nbinom2,data=df.all.larv.long)
summary(larv.mod.nb2)
##  Family: nbinom2  ( log )
## Formula:          
## value ~ offset(log(sum)) + Microbe_treatment * Food_type + (1 |  
##     variable) + (1 | variable:Microbe_treatment) + (1 | variable:Food_type) +  
##     (1 | variable:Microbe_treatment:Food_type)
## Dispersion:             ~variable + Food_type + Microbe_treatment
## Data: df.all.larv.long
## 
##      AIC      BIC   logLik deviance df.resid 
##  15247.7  15678.5  -7550.9  15101.7     2627 
## 
## Random effects:
## 
## Conditional model:
##  Groups                               Name        Variance  Std.Dev.
##  variable                             (Intercept) 1.099e-06 0.001048
##  variable:Microbe_treatment           (Intercept) 1.389e+01 3.727191
##  variable:Food_type                   (Intercept) 2.584e+01 5.082913
##  variable:Microbe_treatment:Food_type (Intercept) 1.727e+00 1.313982
## Number of obs: 2700, groups:  
## variable, 54; variable:Microbe_treatment, 270; variable:Food_type, 108; variable:Microbe_treatment:Food_type, 540
## 
## Conditional model:
##                                           Estimate Std. Error z value Pr(>|z|)
## (Intercept)                                -5.2123     0.9771  -5.334 9.58e-08
## Microbe_treatmentECO                      -10.3785     1.1809  -8.789  < 2e-16
## Microbe_treatmentEMO                       -3.4299     0.9398  -3.649 0.000263
## Microbe_treatmentMEM                       -1.0148     0.9043  -1.122 0.261784
## Microbe_treatmentMMO                       -5.4497     1.0973  -4.967 6.81e-07
## Food_typeLarval food                       -5.6533     1.1879  -4.759 1.95e-06
## Microbe_treatmentECO:Food_typeLarval food  -1.2564     1.3606  -0.923 0.355788
## Microbe_treatmentEMO:Food_typeLarval food   4.3010     0.8108   5.305 1.13e-07
## Microbe_treatmentMEM:Food_typeLarval food   1.6533     0.7612   2.172 0.029853
## Microbe_treatmentMMO:Food_typeLarval food   1.8631     0.9013   2.067 0.038718
##                                              
## (Intercept)                               ***
## Microbe_treatmentECO                      ***
## Microbe_treatmentEMO                      ***
## Microbe_treatmentMEM                         
## Microbe_treatmentMMO                      ***
## Food_typeLarval food                      ***
## Microbe_treatmentECO:Food_typeLarval food    
## Microbe_treatmentEMO:Food_typeLarval food ***
## Microbe_treatmentMEM:Food_typeLarval food *  
## Microbe_treatmentMMO:Food_typeLarval food *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Dispersion model:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)          -3.243182   0.350061  -9.265  < 2e-16 ***
## variableASV0002      -0.591657   0.452883  -1.306 0.191408    
## variableASV0003       1.985713   0.401882   4.941 7.77e-07 ***
## variableASV0004      -1.074704   0.554580  -1.938 0.052639 .  
## variableASV0005       0.066652   0.508792   0.131 0.895775    
## variableASV0006       2.174526   0.438935   4.954 7.27e-07 ***
## variableASV0007       1.385602   0.440130   3.148 0.001643 ** 
## variableASV0009       0.969397   0.432365   2.242 0.024956 *  
## variableASV0010       0.007023   0.450935   0.016 0.987574    
## variableASV0011       1.343682   0.501763   2.678 0.007408 ** 
## variableASV0012      -0.493686   0.508836  -0.970 0.331934    
## variableASV0013       4.735712   0.524659   9.026  < 2e-16 ***
## variableASV0016      -0.193298   0.515311  -0.375 0.707579    
## variableASV0019       0.494949   0.496496   0.997 0.318820    
## variableASV0024       1.352986   0.436170   3.102 0.001922 ** 
## variableASV0027       2.222017   0.485900   4.573 4.81e-06 ***
## variableASV0030      -2.229292   0.576144  -3.869 0.000109 ***
## variableASV0032       4.293441   0.517095   8.303  < 2e-16 ***
## variableASV0033       0.846227   0.411874   2.055 0.039920 *  
## variableASV0043       0.310331   0.567328   0.547 0.584375    
## variableASV0053       1.535773   0.463201   3.316 0.000915 ***
## variableASV0054      -0.727097   0.478087  -1.521 0.128298    
## variableASV0055       2.429137   0.489142   4.966 6.83e-07 ***
## variableASV0058       2.095440   0.500062   4.190 2.79e-05 ***
## variableASV0059       1.057266   0.450351   2.348 0.018892 *  
## variableASV0060       0.889673   0.489455   1.818 0.069113 .  
## variableASV0061       0.565403   0.438701   1.289 0.197463    
## variableASV0066       0.250592   0.500883   0.500 0.616864    
## variableASV0070       2.241231   0.536312   4.179 2.93e-05 ***
## variableASV0079       2.149780   0.425556   5.052 4.38e-07 ***
## variableASV0083       1.649553   0.423276   3.897 9.73e-05 ***
## variableASV0085      -0.793906   0.578679  -1.372 0.170086    
## variableASV0087       3.098546   0.494511   6.266 3.71e-10 ***
## variableASV0099      -0.147528   0.446542  -0.330 0.741114    
## variableASV0100       1.045376   0.506538   2.064 0.039040 *  
## variableASV0102      -0.489850   0.572942  -0.855 0.392565    
## variableASV0109       0.651864   0.620928   1.050 0.293799    
## variableASV0125       1.137772   0.495724   2.295 0.021723 *  
## variableASV0130       0.213252   0.445147   0.479 0.631897    
## variableASV0134      -0.683177   0.520415  -1.313 0.189266    
## variableASV0138       1.987035   0.525266   3.783 0.000155 ***
## variableASV0149       3.648389   0.495782   7.359 1.85e-13 ***
## variableASV0150       1.336130   0.512202   2.609 0.009091 ** 
## variableASV0152       0.804509   0.569131   1.414 0.157486    
## variableASV0164       0.887599   0.477293   1.860 0.062934 .  
## variableASV0210      -1.296908   0.576441  -2.250 0.024458 *  
## variableASV0212       0.758784   0.591754   1.282 0.199751    
## variableASV0218       0.115050   0.581424   0.198 0.843141    
## variableASV0220       0.392096   0.625029   0.627 0.530446    
## variableASV0225       0.452762   0.647796   0.699 0.484598    
## variableASV0262       0.274559   0.787801   0.349 0.727455    
## variableASV0276       0.343902   0.625959   0.549 0.582730    
## variableASV0289       1.540508   0.559074   2.755 0.005861 ** 
## variableASV0311       1.767208   0.540069   3.272 0.001067 ** 
## Food_typeLarval food  0.967602   0.130601   7.409 1.27e-13 ***
## Microbe_treatmentECO  2.851151   0.397070   7.180 6.95e-13 ***
## Microbe_treatmentEMO  0.355197   0.147717   2.405 0.016191 *  
## Microbe_treatmentMEM  0.523868   0.136667   3.833 0.000127 ***
## Microbe_treatmentMMO -0.143964   0.205901  -0.699 0.484434    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
AICtab(larv.mod.nb1,larv.mod.nb2,larv.mod.bin)
##              dAIC  df
## larv.mod.nb1   0.0 73
## larv.mod.nb2 571.6 73
## larv.mod.bin 907.0 15
larv.mod.bin.resid <- simulateResiduals(fittedModel = larv.mod.bin, plot = T)

larv.mod.nb1.resid <- simulateResiduals(fittedModel = larv.mod.nb1, plot = T)

larv.mod.nb2.resid <- simulateResiduals(fittedModel = larv.mod.nb2, plot = T)

##dispersal formula checks
larv.mod.nb1.disnovar <- glmmTMB(value~offset(log(sum))+Microbe_treatment*Food_type+(1|variable)+(1|variable:Microbe_treatment)+(1|variable:Food_type)+(1|variable:Microbe_treatment:Food_type),family=nbinom1,dispformula=~Microbe_treatment+Food_type,data=df.all.larv.long)

anova(larv.mod.nb1.disnovar,larv.mod.nb1) #sig***
## Data: df.all.larv.long
## Models:
## larv.mod.nb1.disnovar: value ~ offset(log(sum)) + Microbe_treatment * Food_type + (1 | , zi=~0, disp=~Microbe_treatment + Food_type
## larv.mod.nb1.disnovar:     variable) + (1 | variable:Microbe_treatment) + (1 | variable:Food_type) + , zi=~0, disp=~variable + Food_type + Microbe_treatment
## larv.mod.nb1.disnovar:     (1 | variable:Microbe_treatment:Food_type), zi=~0, disp=~Microbe_treatment + Food_type
## larv.mod.nb1: value ~ offset(log(sum)) + Microbe_treatment * Food_type + (1 | , zi=~0, disp=~variable + Food_type + Microbe_treatment
## larv.mod.nb1:     variable) + (1 | variable:Microbe_treatment) + (1 | variable:Food_type) + , zi=~0, disp=~Microbe_treatment + Food_type
## larv.mod.nb1:     (1 | variable:Microbe_treatment:Food_type), zi=~0, disp=~variable + Food_type + Microbe_treatment
##                       Df   AIC   BIC  logLik deviance Chisq Chi Df Pr(>Chisq)
## larv.mod.nb1.disnovar 20 14775 14893 -7367.4    14735                        
## larv.mod.nb1          73 14676 15107 -7265.0    14530 204.7     53  < 2.2e-16
##                          
## larv.mod.nb1.disnovar    
## larv.mod.nb1          ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##removing 3 way interaction
larv.mod.nb1.nofm <- update(larv.mod.nb1,.~.-(1|variable:Microbe_treatment:Food_type))
summary(larv.mod.nb1.nofm)
##  Family: nbinom1  ( log )
## Formula:          
## value ~ Microbe_treatment + Food_type + (1 | variable) + (1 |  
##     variable:Microbe_treatment) + (1 | variable:Food_type) +  
##     Microbe_treatment:Food_type + offset(log(sum))
## Dispersion:             ~variable + Food_type + Microbe_treatment
## Data: df.all.larv.long
## 
##      AIC      BIC   logLik deviance df.resid 
##  14678.5  15103.4  -7267.3  14534.5     2628 
## 
## Random effects:
## 
## Conditional model:
##  Groups                     Name        Variance  Std.Dev. 
##  variable                   (Intercept) 1.312e-07 0.0003623
##  variable:Microbe_treatment (Intercept) 4.029e+00 2.0072916
##  variable:Food_type         (Intercept) 6.676e+00 2.5837467
## Number of obs: 2700, groups:  
## variable, 54; variable:Microbe_treatment, 270; variable:Food_type, 108
## 
## Conditional model:
##                                           Estimate Std. Error z value Pr(>|z|)
## (Intercept)                                -6.4678     0.4802 -13.469  < 2e-16
## Microbe_treatmentECO                       -2.8162     0.9564  -2.945 0.003233
## Microbe_treatmentEMO                       -0.5032     0.4541  -1.108 0.267817
## Microbe_treatmentMEM                       -0.1419     0.4287  -0.331 0.740638
## Microbe_treatmentMMO                       -2.2053     0.4972  -4.435 9.19e-06
## Food_typeLarval food                       -1.2519     0.5825  -2.149 0.031613
## Microbe_treatmentECO:Food_typeLarval food  -2.6340     0.5224  -5.042 4.61e-07
## Microbe_treatmentEMO:Food_typeLarval food   0.9767     0.2841   3.438 0.000586
## Microbe_treatmentMEM:Food_typeLarval food   0.7184     0.2366   3.036 0.002395
## Microbe_treatmentMMO:Food_typeLarval food   1.1724     0.2773   4.228 2.36e-05
##                                              
## (Intercept)                               ***
## Microbe_treatmentECO                      ** 
## Microbe_treatmentEMO                         
## Microbe_treatmentMEM                         
## Microbe_treatmentMMO                      ***
## Food_typeLarval food                      *  
## Microbe_treatmentECO:Food_typeLarval food ***
## Microbe_treatmentEMO:Food_typeLarval food ***
## Microbe_treatmentMEM:Food_typeLarval food ** 
## Microbe_treatmentMMO:Food_typeLarval food ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Dispersion model:
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)           6.85166    0.46809  14.638  < 2e-16 ***
## variableASV0002       0.04264    0.85391   0.050 0.960178    
## variableASV0003       1.17327    0.53168   2.207 0.027334 *  
## variableASV0004       1.15220    0.98560   1.169 0.242390    
## variableASV0005       0.38587    0.72643   0.531 0.595292    
## variableASV0006       1.35183    0.56042   2.412 0.015857 *  
## variableASV0007       0.02670    0.57158   0.047 0.962747    
## variableASV0009       0.61519    0.57860   1.063 0.287675    
## variableASV0010      -0.07145    0.63805  -0.112 0.910843    
## variableASV0011      -0.92455    0.62355  -1.483 0.138146    
## variableASV0012       0.20960    0.81152   0.258 0.796190    
## variableASV0013       1.14764    0.65803   1.744 0.081149 .  
## variableASV0016       0.52626    0.76301   0.690 0.490370    
## variableASV0019       1.58854    0.66730   2.381 0.017287 *  
## variableASV0024      -0.33074    0.58623  -0.564 0.572629    
## variableASV0027       0.36266    0.63071   0.575 0.565290    
## variableASV0030      -0.14127    1.08752  -0.130 0.896644    
## variableASV0032       0.21191    0.70540   0.300 0.763861    
## variableASV0033       1.08423    0.60075   1.805 0.071105 .  
## variableASV0043       3.09257    0.84087   3.678 0.000235 ***
## variableASV0053       0.10922    0.61851   0.177 0.859838    
## variableASV0054       0.26393    0.83115   0.318 0.750833    
## variableASV0055       1.21154    0.68891   1.759 0.078641 .  
## variableASV0058       1.41797    0.70361   2.015 0.043875 *  
## variableASV0059      -0.71776    0.60315  -1.190 0.234036    
## variableASV0060       1.82663    0.75287   2.426 0.015256 *  
## variableASV0061       1.03000    0.64992   1.585 0.113010    
## variableASV0066      -1.42025    0.68413  -2.076 0.037896 *  
## variableASV0070      -0.27711    0.81396  -0.340 0.733523    
## variableASV0079      -0.95055    0.57573  -1.651 0.098736 .  
## variableASV0083      -0.77493    0.58019  -1.336 0.181665    
## variableASV0085      -0.18678    0.95542  -0.195 0.845007    
## variableASV0087       0.27474    0.64468   0.426 0.669991    
## variableASV0099       0.24630    0.72144   0.341 0.732797    
## variableASV0100      -0.47850    0.65498  -0.731 0.465046    
## variableASV0102      -1.57001    0.85924  -1.827 0.067667 .  
## variableASV0109       0.32935    1.00606   0.327 0.743391    
## variableASV0125       1.08067    0.76783   1.407 0.159302    
## variableASV0130      -0.01441    0.68679  -0.021 0.983255    
## variableASV0134       0.73435    0.90294   0.813 0.416053    
## variableASV0138      -1.13267    0.73269  -1.546 0.122123    
## variableASV0149      -2.29236    0.61342  -3.737 0.000186 ***
## variableASV0150      -1.53318    0.68541  -2.237 0.025293 *  
## variableASV0152       0.61808    0.89867   0.688 0.491598    
## variableASV0164      -1.33252    0.65650  -2.030 0.042382 *  
## variableASV0210       0.17882    1.06693   0.168 0.866893    
## variableASV0212       0.35378    0.95374   0.371 0.710678    
## variableASV0218      -1.30772    0.84022  -1.556 0.119611    
## variableASV0220       0.81224    1.04553   0.777 0.437237    
## variableASV0225       0.02016    1.06207   0.019 0.984859    
## variableASV0262      -0.56844    1.04866  -0.542 0.587775    
## variableASV0276       0.43079    1.02047   0.422 0.672915    
## variableASV0289      -1.44067    0.79490  -1.812 0.069925 .  
## variableASV0311      -1.22539    0.76813  -1.595 0.110647    
## Food_typeLarval food  0.96112    0.21909   4.387 1.15e-05 ***
## Microbe_treatmentECO -1.07382    0.82793  -1.297 0.194634    
## Microbe_treatmentEMO -0.07345    0.22125  -0.332 0.739912    
## Microbe_treatmentMEM -0.61978    0.19197  -3.229 0.001244 ** 
## Microbe_treatmentMMO  0.21159    0.26152   0.809 0.418463    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(larv.mod.nb1,larv.mod.nb1.nofm) #sig **
## Data: df.all.larv.long
## Models:
## larv.mod.nb1.nofm: value ~ Microbe_treatment + Food_type + (1 | variable) + (1 | , zi=~0, disp=~variable + Food_type + Microbe_treatment
## larv.mod.nb1.nofm:     variable:Microbe_treatment) + (1 | variable:Food_type) + , zi=~0, disp=~variable + Food_type + Microbe_treatment
## larv.mod.nb1.nofm:     Microbe_treatment:Food_type + offset(log(sum)), zi=~0, disp=~variable + Food_type + Microbe_treatment
## larv.mod.nb1: value ~ offset(log(sum)) + Microbe_treatment * Food_type + (1 | , zi=~0, disp=~variable + Food_type + Microbe_treatment
## larv.mod.nb1:     variable) + (1 | variable:Microbe_treatment) + (1 | variable:Food_type) + , zi=~0, disp=~variable + Food_type + Microbe_treatment
## larv.mod.nb1:     (1 | variable:Microbe_treatment:Food_type), zi=~0, disp=~variable + Food_type + Microbe_treatment
##                   Df   AIC   BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)  
## larv.mod.nb1.nofm 72 14678 15103 -7267.3    14534                           
## larv.mod.nb1      73 14676 15107 -7265.0    14530 4.4686      1    0.03452 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# ##removing just food
# larv.mod.nb1.nofood <- update(larv.mod.nb1,.~.-(1|variable:Food_type))
# anova(larv.mod.nb1,larv.mod.nb1.nofood) #super sig***
# 
# ##removing just microbes
# larv.mod.nb1.nomicr <- update(larv.mod.nb1,.~.-(1|variable:Microbe_treatment))
# anova(larv.mod.nb1,larv.mod.nb1.nomicr) #super sig***

##checking out nbinom2
larv.mod.nb2.nofm <- update(larv.mod.nb2,.~.-(1|variable:Microbe_treatment:Food_type))
anova(larv.mod.nb2,larv.mod.nb2.nofm)
## Data: df.all.larv.long
## Models:
## larv.mod.nb2.nofm: value ~ Microbe_treatment + Food_type + (1 | variable) + (1 | , zi=~0, disp=~variable + Food_type + Microbe_treatment
## larv.mod.nb2.nofm:     variable:Microbe_treatment) + (1 | variable:Food_type) + , zi=~0, disp=~variable + Food_type + Microbe_treatment
## larv.mod.nb2.nofm:     Microbe_treatment:Food_type + offset(log(sum)), zi=~0, disp=~variable + Food_type + Microbe_treatment
## larv.mod.nb2: value ~ offset(log(sum)) + Microbe_treatment * Food_type + (1 | , zi=~0, disp=~variable + Food_type + Microbe_treatment
## larv.mod.nb2:     variable) + (1 | variable:Microbe_treatment) + (1 | variable:Food_type) + , zi=~0, disp=~variable + Food_type + Microbe_treatment
## larv.mod.nb2:     (1 | variable:Microbe_treatment:Food_type), zi=~0, disp=~variable + Food_type + Microbe_treatment
##                   Df   AIC   BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)   
## larv.mod.nb2.nofm 72 15256 15681 -7556.1    15112                            
## larv.mod.nb2      73 15248 15678 -7550.9    15102 10.524      1   0.001179 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##sig**

##checking out binomial
larv.mod.bin.nofm <- update(larv.mod.bin,.~.-(1|variable:Microbe_treatment:Food_type))
summary(larv.mod.bin.nofm)
##  Family: binomial  ( logit )
## Formula:          
## cbind(value, sum - value) ~ Microbe_treatment + Food_type + (1 |  
##     variable) + (1 | variable:Microbe_treatment) + (1 | variable:Food_type) +  
##     (1 | rowid) + Microbe_treatment:Food_type
## Data: df.all.larv.long
## 
##      AIC      BIC   logLik deviance df.resid 
##  15601.7  15684.3  -7786.9  15573.7     2686 
## 
## Random effects:
## 
## Conditional model:
##  Groups                     Name        Variance  Std.Dev.
##  variable                   (Intercept) 1.307e-06 0.001143
##  variable:Microbe_treatment (Intercept) 3.089e+01 5.557434
##  variable:Food_type         (Intercept) 4.042e+01 6.357299
##  rowid                      (Intercept) 1.406e+01 3.749305
## Number of obs: 2700, groups:  
## variable, 54; variable:Microbe_treatment, 270; variable:Food_type, 108; rowid, 2700
## 
## Conditional model:
##                                           Estimate Std. Error z value Pr(>|z|)
## (Intercept)                               -11.3323     1.2193  -9.294  < 2e-16
## Microbe_treatmentECO                       -9.0397     1.3069  -6.917 4.62e-12
## Microbe_treatmentEMO                       -2.6017     1.2026  -2.163  0.03051
## Microbe_treatmentMEM                        0.4579     1.1676   0.392  0.69494
## Microbe_treatmentMMO                       -8.3851     1.1753  -7.135 9.70e-13
## Food_typeLarval food                       -5.8334     1.4727  -3.961 7.46e-05
## Microbe_treatmentECO:Food_typeLarval food  -2.8321     1.4541  -1.948  0.05146
## Microbe_treatmentEMO:Food_typeLarval food   3.5878     0.7640   4.696 2.65e-06
## Microbe_treatmentMEM:Food_typeLarval food   1.8169     0.7501   2.422  0.01543
## Microbe_treatmentMMO:Food_typeLarval food   2.9328     0.8917   3.289  0.00101
##                                              
## (Intercept)                               ***
## Microbe_treatmentECO                      ***
## Microbe_treatmentEMO                      *  
## Microbe_treatmentMEM                         
## Microbe_treatmentMMO                      ***
## Food_typeLarval food                      ***
## Microbe_treatmentECO:Food_typeLarval food .  
## Microbe_treatmentEMO:Food_typeLarval food ***
## Microbe_treatmentMEM:Food_typeLarval food *  
## Microbe_treatmentMMO:Food_typeLarval food ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(larv.mod.bin,larv.mod.bin.nofm) #sig***
## Data: df.all.larv.long
## Models:
## larv.mod.bin.nofm: cbind(value, sum - value) ~ Microbe_treatment + Food_type + (1 | , zi=~0, disp=~1
## larv.mod.bin.nofm:     variable) + (1 | variable:Microbe_treatment) + (1 | variable:Food_type) + , zi=~0, disp=~1
## larv.mod.bin.nofm:     (1 | rowid) + Microbe_treatment:Food_type, zi=~0, disp=~1
## larv.mod.bin: cbind(value, sum - value) ~ Microbe_treatment * Food_type + (1 | , zi=~0, disp=~1
## larv.mod.bin:     variable) + (1 | variable:Microbe_treatment) + (1 | variable:Food_type) + , zi=~0, disp=~1
## larv.mod.bin:     (1 | variable:Microbe_treatment:Food_type) + (1 | rowid), zi=~0, disp=~1
##                   Df   AIC   BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
## larv.mod.bin.nofm 14 15602 15684 -7786.9    15574                             
## larv.mod.bin      15 15583 15672 -7776.6    15553 20.617      1   5.61e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
###save conditional modes
larv.mod.ranef <- as.data.frame(ranef(larv.mod.nb1))
write.csv(larv.mod.ranef, "larv.mod.ranef.all.csv")
saveRDS(larv.mod.nb1,file="larv.mod.nb1.all.rds")

# library("ggstats")
# ggcoef_multicomponents(larv.mod.nb1)

3.4.1 Nbinom1 vs. nbinom2 things

##note: some errors/warnings if the treatments aren't factors for some reason
mvtab <- ddply(df.all.larv.long,
               .(variable:Microbe_treatment:Food_type),
               summarise,
               callmean=mean(value),
               callvar=var(value))

q1 <- qplot(callmean,callvar,data=mvtab)
print(q1+
        ## linear (quasi-Poisson/NB1) fit
        geom_smooth(method="lm",formula=y~x-1,colour="pink")+
        ## smooth (loess)
        geom_smooth(colour="red")+
        ## semi-quadratic (NB2/LNP)
        geom_smooth(method="lm",formula=y~I(x^2)+offset(x)-1,colour="purple")+
        ## Poisson (v=m)
        geom_abline(a=0,b=1,lty=2))#+
## Warning in geom_abline(a = 0, b = 1, lty = 2): Ignoring unknown parameters: `a`
## and `b`
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

        #ylim(0,1.75e8)+
        #xlim(0,20000)